diff --git a/pygrt/C_extension/include/grt/common/RT_matrix.h b/pygrt/C_extension/include/grt/common/RT_matrix.h index d92b01c..ed6478e 100644 --- a/pygrt/C_extension/include/grt/common/RT_matrix.h +++ b/pygrt/C_extension/include/grt/common/RT_matrix.h @@ -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 diff --git a/pygrt/C_extension/include/grt/common/const.h b/pygrt/C_extension/include/grt/common/const.h index def2bd5..a645846 100755 --- a/pygrt/C_extension/include/grt/common/const.h +++ b/pygrt/C_extension/include/grt/common/const.h @@ -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 { diff --git a/pygrt/C_extension/include/grt/common/model.h b/pygrt/C_extension/include/grt/common/model.h index 6dbfe67..416ad0a 100755 --- a/pygrt/C_extension/include/grt/common/model.h +++ b/pygrt/C_extension/include/grt/common/model.h @@ -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]; @@ -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); /** * 从模型文件中判断各个量的大致精度（字符串长度），以确定浮点数输出位数 diff --git a/pygrt/C_extension/include/grt/dynamic/layer.h b/pygrt/C_extension/include/grt/dynamic/layer.h index ca4e938..75af04b 100755 --- a/pygrt/C_extension/include/grt/dynamic/layer.h +++ b/pygrt/C_extension/include/grt/dynamic/layer.h @@ -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); /** diff --git a/pygrt/C_extension/include/grt/static/static_layer.h b/pygrt/C_extension/include/grt/static/static_layer.h index 9115697..7be2e83 100644 --- a/pygrt/C_extension/include/grt/static/static_layer.h +++ b/pygrt/C_extension/include/grt/static/static_layer.h @@ -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) diff --git a/pygrt/C_extension/src/common/RT_matrix.c b/pygrt/C_extension/src/common/RT_matrix.c index 9be9896..d0f7d23 100644 --- a/pygrt/C_extension/src/common/RT_matrix.c +++ b/pygrt/C_extension/src/common/RT_matrix.c @@ -13,12 +13,31 @@ #include #include #include +#include #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); diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index 853e0ca..64378b3 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -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; diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index c20e23f..2ccce70 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -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; @@ -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); @@ -289,7 +298,7 @@ printf("\n" " -R,[,...] -O [-H/] \n" " [-L] [-C[d|p|n]] [-E[/]] \n" " [-K[+k][+s][+e][+v]]\n" -" [-P] [-Ge|v|h|s] \n" +" [-P] [-Ge|v|h|s] [-Bf|F|r|R|h|H]\n" " [-S[,,...]] [-e] [-s]\n" "\n\n" "Options:\n" @@ -397,6 +406,13 @@ printf("\n" " : Horizontal Force (HF)\n" " : 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[,,...]\n" " Frequency (index) of statsfile in wavenumber\n" " integration to be output, require 0 <= i <= nf-1,\n" @@ -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; @@ -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 @@ -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][+n][+a] case 'N': Ctrl->N.active = true; @@ -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){ diff --git a/pygrt/C_extension/src/dynamic/layer.c b/pygrt/C_extension/src/dynamic/layer.c index b34cf77..7e31703 100644 --- a/pygrt/C_extension/src/dynamic/layer.c +++ b/pygrt/C_extension/src/dynamic/layer.c @@ -29,37 +29,118 @@ T N##2 = mod1d->N[iy];\ -void grt_topfree_RU(GRT_MODEL1D *mod1d) +static int __freeBound_R(cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML) { - cplx_t cbcb0 = mod1d->cbcb[0]; if(cbcb0 != 0.0){ // 固体表面 // 公式(5.3.10-14) cplx_t Delta = 0.0; cplx_t cbcb02 = 0.25*cbcb0*cbcb0; - cplx_t xa0 = mod1d->xa[0]; - cplx_t xb0 = mod1d->xb[0]; - // 对公式(5.3.10-14)进行重新整理，对浮点数友好一些 Delta = -1.0 + xa0*xb0 + cbcb0 - cbcb02; if(Delta == 0.0){ - mod1d->M_top.stats = GRT_INVERSE_FAILURE; - return; + return GRT_INVERSE_FAILURE; } Delta = 1.0 / Delta; - mod1d->M_top.RU[0][0] = (1.0 + xa0*xb0 - cbcb0 + cbcb02) * Delta; - mod1d->M_top.RU[0][1] = 2.0 * xb0 * (1.0 - 0.5*cbcb0) * Delta; - mod1d->M_top.RU[1][0] = 2.0 * xa0 * (1.0 - 0.5*cbcb0) * Delta; - mod1d->M_top.RU[1][1] = mod1d->M_top.RU[0][0]; - mod1d->M_top.RUL = 1.0; + M[0][0] = (1.0 + xa0*xb0 - cbcb0 + cbcb02) * Delta; + M[0][1] = 2.0 * xb0 * (1.0 - 0.5*cbcb0) * Delta * adsgn; + M[1][0] = 2.0 * xa0 * (1.0 - 0.5*cbcb0) * Delta * adsgn; + M[1][1] = M[0][0]; + *ML = 1.0; } else { // 液体表面 - mod1d->M_top.RU[0][0] = -1.0; - mod1d->M_top.RU[1][1] = mod1d->M_top.RU[0][1] = mod1d->M_top.RU[1][0] = 0.0; + M[0][0] = -1.0; + M[1][1] = M[0][1] = M[1][0] = 0.0; + *ML = 0.0; + } + return GRT_INVERSE_SUCCESS; +} + +static int __rigidBound_R(cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML) +{ + if(cbcb0 != 0.0){ + // 固体表面 + cplx_t Delta = 0.0; + + Delta = xa0*xb0 - 1.0; + if(Delta == 0.0){ + return GRT_INVERSE_FAILURE; + } + Delta = 1.0 / Delta; + M[0][0] = (xa0*xb0 + 1.0) * Delta; + M[0][1] = 2.0 * xb0 * Delta * adsgn; + M[1][0] = 2.0 * xa0 * Delta * adsgn; + M[1][1] = M[0][0]; + *ML = - 1.0; + } + else { + // 液体表面 + M[0][0] = 1.0; + M[1][1] = M[0][1] = M[1][0] = 0.0; + *ML = 0.0; + } + return GRT_INVERSE_SUCCESS; +} + +void grt_topbound_RU(GRT_MODEL1D *mod1d) +{ + cplx_t cbcb = mod1d->cbcb[0]; + cplx_t xa = mod1d->xa[0]; + cplx_t xb = mod1d->xb[0]; + if(mod1d->topbound == GRT_BOUND_FREE){ + mod1d->M_top.stats = __freeBound_R(cbcb, xa, xb, 1.0, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_RIGID){ + mod1d->M_top.stats = __rigidBound_R(cbcb, xa, xb, 1.0, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_HALFSPACE){ + memset(mod1d->M_top.RU, 0, sizeof(cplx_t)*4); mod1d->M_top.RUL = 0.0; + mod1d->M_top.stats = GRT_INVERSE_SUCCESS; } + else{ + GRTRaiseError("Wrong execution."); + } + // RU 不需要时延 +} + +void grt_botbound_RD(GRT_MODEL1D *mod1d) +{ + size_t nlay = mod1d->n; + cplx_t cbcb = mod1d->cbcb[nlay-2]; + cplx_t xa = mod1d->xa[nlay-2]; + cplx_t xb = mod1d->xb[nlay-2]; + if(mod1d->botbound == GRT_BOUND_FREE){ + mod1d->M_bot.stats = __freeBound_R(cbcb, xa, xb, -1.0, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_RIGID){ + mod1d->M_bot.stats = __rigidBound_R(cbcb, xa, xb, -1.0, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_HALFSPACE){ + grt_RT_matrix_PSV(mod1d, nlay-1, &mod1d->M_bot); + grt_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot); + } + else{ + GRTRaiseError("Wrong execution."); + } + + // 时延 RD + real_t k = mod1d->k; + real_t thk = mod1d->Thk[nlay-2]; + + cplx_t exa, exb, ex2a, ex2b, exab; + exa = exp(- k*thk*xa); + exb = exp(- k*thk*xb); + ex2a = exa * exa; + ex2b = exb * exb; + exab = exa * exb; + + mod1d->M_bot.RD[0][0] *= ex2a; mod1d->M_bot.RD[0][1] *= exab; + mod1d->M_bot.RD[1][0] *= exab; mod1d->M_bot.RD[1][1] *= ex2b; + + mod1d->M_bot.RDL *= ex2b; } diff --git a/pygrt/C_extension/src/integral/kernel_template.c_ b/pygrt/C_extension/src/integral/kernel_template.c_ index caca5de..1448a52 100644 --- a/pygrt/C_extension/src/integral/kernel_template.c_ +++ b/pygrt/C_extension/src/integral/kernel_template.c_ @@ -47,7 +47,8 @@ #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_topbound_RU CONCAT_PREFIX(topbound_RU) +#define __grt_botbound_RD CONCAT_PREFIX(botbound_RD) #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) @@ -72,14 +73,23 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) RT_MATRIX *M_BL = &mod1d->M_BL; RT_MATRIX *M_RS = &mod1d->M_RS; RT_MATRIX *M_FA = &mod1d->M_FA; - // 自由表面 + // 顶底界面 RT_MATRIX *M_top = &mod1d->M_top; + RT_MATRIX *M_bot = &mod1d->M_bot; + + grt_reset_RT_matrix(M_BL); + grt_reset_RT_matrix(M_RS); + grt_reset_RT_matrix(M_FA); + grt_reset_RT_matrix(M_top); + grt_reset_RT_matrix(M_bot); // 从顶到底进行矩阵递推, 公式(5.5.3) - for(size_t iy=1; iyn; ++iy){ // 因为n>=3, 故一定会进入该循环 - + size_t maxidx = (imax+1 == mod1d->n)? mod1d->n : mod1d->n - 1; + for(size_t iy = 1; iy < maxidx; ++iy) // 因为n>=3, 故一定会进入该循环 + { // 定义物理层内的反射透射系数矩阵 - grt_init_RT_matrix(M); + RT_MATRIX *M = &(RT_MATRIX){}; + grt_reset_RT_matrix(M); // 只有动态解才可以使用这个 if 简化， // 对于静态解，即使是震源层、接收层这种虚拟层位也需要显式计算R/T矩阵 @@ -143,11 +153,20 @@ void __GRT_MATRIX_FUNC__(GRT_MODEL1D *mod1d, const real_t k) //=================================================================================== // 递推RU_FA - __grt_topfree_RU(mod1d); + __grt_topbound_RU(mod1d); if(M_top->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; grt_recursion_RU(M_top, M_FA, M_FA); if(M_FA->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + // 递推RD_BL, 在模型定义阶段已检查震源/台站的虚拟层成为底层时，只能是半无限空间 + if(imax+1 < mod1d->n){ // 此 if 仅减少不必要的计算 + __grt_botbound_RD(mod1d); + if(M_bot->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + grt_recursion_RD(M_BL, M_bot, M_BL); + if(M_BL->stats==GRT_INVERSE_FAILURE) goto BEFORE_RETURN; + } + + BEFORE_RETURN: return; } diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index 29a92f9..4b6270d 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -41,6 +41,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; @@ -143,9 +149,9 @@ printf("\n" "\n\n" "Usage:\n" "----------------------------------------------------------------\n" -" grt static greenfn -M -D/ -X// \n" +" grt static greenfn -M -D/ -X// \n" " -Y// -O [-L] [-C[d|p|n]] \n" -" [-K[+k][+e]] [-S] [-e]\n" +" [-Bf|F|r|R|h|H] [-K[+k][+e]] [-S] [-e]\n" "\n\n" "Options:\n" "----------------------------------------------------------------\n" @@ -198,6 +204,13 @@ printf("\n" " + n: None.\n" " Default use +cd when fabs(depsrc-deprcv) <= %.1f.\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf( "\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" " -K[+k][+e]\n" " Several parameters designed to define the\n" " behavior in wavenumber integration. The upper\n" @@ -236,11 +249,14 @@ 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->K.k0 = GRT_GREENFN_K_K0; int opt; - while ((opt = getopt(argc, argv, ":M:D:L:C:K:X:Y:O:Seh")) != -1) { + while ((opt = getopt(argc, argv, ":M:D:B:L:C:K:X:Y:O:Seh")) != -1) { switch (opt) { // 模型路径，其中每行分别为 // 厚度(km) Vp(km/s) Vs(km/s) Rho(g/cm^3) Qp Qs @@ -270,6 +286,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; + // 波数积分间隔 -L[][+l][+a][+o] case 'L': Ctrl->L.active = true; @@ -482,6 +516,9 @@ int static_greenfn_main(int argc, char **argv){ } GRT_MODEL1D *mod1d = Ctrl->M.mod1d; + // 边界条件 + grt_set_mod1d_boundary(mod1d, Ctrl->B.topbound, Ctrl->B.botbound); + // 判断是否要自动使用收敛方法 if( ! Ctrl->C.active && fabs(Ctrl->D.deprcv - Ctrl->D.depsrc) <= GRT_MIN_DEPTH_GAP_SRC_RCV) { Ctrl->C.applyDCM = true; diff --git a/pygrt/C_extension/src/static/static_layer.c b/pygrt/C_extension/src/static/static_layer.c index 687f927..a614a62 100644 --- a/pygrt/C_extension/src/static/static_layer.c +++ b/pygrt/C_extension/src/static/static_layer.c @@ -14,6 +14,7 @@ #include #include #include +#include #include "grt/static/static_layer.h" #include "grt/common/model.h" @@ -25,14 +26,87 @@ T N##2 = mod1d->N[iy];\ -void grt_static_topfree_RU(GRT_MODEL1D *mod1d) +static int __freeBound_R(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2], cplx_t *ML) { - cplx_t delta1 = mod1d->delta[0]; - // 公式(6.3.12) - mod1d->M_top.RU[0][0] = mod1d->M_top.RU[1][1] = 0.0; - mod1d->M_top.RU[0][1] = -delta1; - mod1d->M_top.RU[1][0] = -1.0/delta1; - mod1d->M_top.RUL = 1.0; + if(d <= 0.0){ + // 公式(6.3.12) + M[0][0] = M[1][1] = 0.0; + M[0][1] = -delta1; + M[1][0] = -1.0/delta1; + } else { + M[0][0] = M[1][1] = -2.0*d*k; + M[0][1] = delta1 * (-4.0*d*d*k*k - 1.0); + M[1][0] = -1.0/delta1; + } + *ML = 1.0; + return GRT_INVERSE_SUCCESS; +} + +static int __rigidBound_R(real_t d, real_t k, cplx_t delta1, cplx_t M[2][2], cplx_t *ML) +{ + if(d <= 0.0){ + // 公式(6.3.12) + M[0][0] = M[1][1] = 0.0; + M[0][1] = M[1][0] = 1.0; + } else { + M[0][0] = M[1][1] = 2.0*delta1*d*k; + M[0][1] = (4.0*GRT_SQUARE(delta1*d*k) + 1.0); + M[1][0] = 1.0; + } + *ML = -1.0; + return GRT_INVERSE_SUCCESS; +} + +void grt_static_topbound_RU(GRT_MODEL1D *mod1d) +{ + cplx_t delta = mod1d->delta[0]; + real_t k = mod1d->k; + if(mod1d->topbound == GRT_BOUND_FREE){ + mod1d->M_top.stats = __freeBound_R(0.0, k, delta, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_RIGID){ + mod1d->M_top.stats = __rigidBound_R(0.0, k, delta, mod1d->M_top.RU, &mod1d->M_top.RUL); + } + else if(mod1d->topbound == GRT_BOUND_HALFSPACE){ + memset(mod1d->M_top.RU, 0, sizeof(cplx_t)*4); + mod1d->M_top.RUL = 0.0; + mod1d->M_top.stats = GRT_INVERSE_SUCCESS; + } + else{ + GRTRaiseError("Wrong execution."); + } + // RU 不需要时延 +} + +void grt_static_botbound_RD(GRT_MODEL1D *mod1d) +{ + size_t nlay = mod1d->n; + cplx_t delta = mod1d->delta[nlay-2]; + real_t thk = mod1d->Thk[nlay-2]; + real_t k = mod1d->k; + if(mod1d->botbound == GRT_BOUND_FREE){ + mod1d->M_bot.stats = __freeBound_R(thk, k, delta, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_RIGID){ + mod1d->M_bot.stats = __rigidBound_R(thk, k, delta, mod1d->M_bot.RD, &mod1d->M_bot.RDL); + } + else if(mod1d->botbound == GRT_BOUND_HALFSPACE){ + grt_static_RT_matrix_PSV(mod1d, nlay-1, &mod1d->M_bot); + grt_static_RT_matrix_SH(mod1d, nlay-1, &mod1d->M_bot); + } + else{ + GRTRaiseError("Wrong execution."); + } + + // 时延 RD + cplx_t ex, ex2; + ex = exp(- k*thk); + ex2 = ex * ex; + + mod1d->M_bot.RD[0][0] *= ex2; mod1d->M_bot.RD[0][1] *= ex2; + mod1d->M_bot.RD[1][0] *= ex2; mod1d->M_bot.RD[1][1] *= ex2; + + mod1d->M_bot.RDL *= ex2; } void grt_static_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d) diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py index 96ae845..62f1ad4 100755 --- a/pygrt/c_interfaces.py +++ b/pygrt/c_interfaces.py @@ -85,6 +85,11 @@ def set_num_threads(n): C_grt_read_mod1d_from_file.restype = POINTER(c_GRT_MODEL1D) C_grt_read_mod1d_from_file.argtypes = [c_char_p, c_double, c_double, c_bool] +C_grt_set_mod1d_boundary = libgrt.grt_set_mod1d_boundary +"""设置模型边界条件并检查底界面""" +C_grt_set_mod1d_boundary.restype = None +C_grt_set_mod1d_boundary.argtypes = [POINTER(c_GRT_MODEL1D), c_int, c_int] + C_grt_free_mod1d = libgrt.grt_free_mod1d """释放C程序中申请的 GRT_MODEL1D 结构体内存""" C_grt_free_mod1d.restype = None diff --git a/pygrt/c_structures.py b/pygrt/c_structures.py index ce17a51..f34aa83 100755 --- a/pygrt/c_structures.py +++ b/pygrt/c_structures.py @@ -123,7 +123,10 @@ class c_GRT_MODEL1D(Structure): ('M_FA', c_RT_MATRIX), ('M_FB', c_RT_MATRIX), + ('topbound', c_int), ('M_top', c_RT_MATRIX), + ('botbound', c_int), + ('M_bot', c_RT_MATRIX), ('R_EV', CPLX * 4), ('R_EVL', CPLX), diff --git a/pygrt/pymod.py b/pygrt/pymod.py index 253008d..37398dc 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -15,7 +15,7 @@ from obspy import read, Stream, Trace, UTCDateTime from scipy.fft import irfft, ifft from obspy.core import AttribDict -from typing import List, Dict, Union +from typing import List, Dict, Union, Literal import tempfile from time import time @@ -33,7 +33,9 @@ class PyModel1D: - def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:bool=False): + def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:bool=False, + topbound:Literal['free', 'rigid', 'halfspace']='free', + botbound:Literal['free', 'rigid', 'halfspace']='halfspace'): ''' Create 1D model instance, and insert the imaginary layer of source and receiver. @@ -41,12 +43,27 @@ def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:b :param depsrc: source depth (km) :param deprcv: receiver depth (km) :param allowLiquid: whether liquid layers are allowed + :param topbound: boundary condition of the top layer + :param botbound: boundary condition of the bottom layer ''' self.depsrc:float = depsrc self.deprcv:float = deprcv self.c_mod1d:c_GRT_MODEL1D self.hasLiquid:bool = allowLiquid # 传入的模型是否有液体层 + self.topbound:str = topbound + self.botbound:str = botbound + + boundDct = { + 'free': 0, + 'rigid': 1, + 'halfspace': 2 + } + + if topbound not in boundDct: + raise ValueError(f"Unsupported topbound={topbound}.") + if botbound not in boundDct: + raise ValueError(f"Unsupported botbound={botbound}.") # 将modarr写入临时数组 with tempfile.NamedTemporaryFile(mode='w', delete=False) as tmpfile: @@ -57,10 +74,15 @@ def __init__(self, modarr0:np.ndarray, depsrc:float, deprcv:float, allowLiquid:b try: c_mod1d_ptr = C_grt_read_mod1d_from_file(tmp_path.encode("utf-8"), depsrc, deprcv, allowLiquid) self.c_mod1d = c_mod1d_ptr.contents # 这部分内存在C中申请，需由C函数释放。占用不多，这里跳过 + C_grt_set_mod1d_boundary(self.c_mod1d, boundDct[topbound], boundDct[botbound]) finally: if os.path.exists(tmp_path): os.unlink(tmp_path) + # 设置边界条件 + self.c_mod1d.topbound = boundDct[topbound] + self.c_mod1d.botbound = boundDct[botbound] + self.isrc = self.c_mod1d.isrc self.ircv = self.c_mod1d.ircv @@ -160,7 +182,6 @@ def _get_grn_spectra( delayT0:float=0.0, delayV0:float=0.0, calc_upar:bool=False, - gf_source=['EX', 'VF', 'HF', 'DC'], statsfile:Union[str,None]=None, statsidxs:Union[np.ndarray,List[int],None]=None, print_runtime:bool=True @@ -439,7 +460,7 @@ def compute_grn( pygrnLst, pygrnLst_uiz, pygrnLst_uir = self._get_grn_spectra( distarr, nt, dt, upsampling_n, freqband, zeta, keepAllFreq, vmin_ref, keps, ampk, k0, Length, filonLength, safilonTol, filonCut, converg_method, - delayT0, delayV0, calc_upar, gf_source, + delayT0, delayV0, calc_upar, statsfile, statsidxs, print_runtime )