From 22ad49846aa7a7276950e2fbd60034de6e37782c Mon Sep 17 00:00:00 2001 From: Dengda98 Date: Thu, 29 Jan 2026 23:05:12 +0800 Subject: [PATCH] REFAC: add a boolean array in struct GRT_MODEL1D to indicate whether each layer is a liquid layer --- pygrt/C_extension/include/grt/common/model.h | 1 + pygrt/C_extension/src/common/model.c | 8 ++- pygrt/C_extension/src/dynamic/grt_greenfn.c | 4 +- pygrt/C_extension/src/dynamic/layer.c | 60 ++++++++++---------- pygrt/c_structures.py | 1 + 5 files changed, 40 insertions(+), 34 deletions(-) diff --git a/pygrt/C_extension/include/grt/common/model.h b/pygrt/C_extension/include/grt/common/model.h index 9042a36..0c73f55 100755 --- a/pygrt/C_extension/include/grt/common/model.h +++ b/pygrt/C_extension/include/grt/common/model.h @@ -38,6 +38,7 @@ typedef struct { real_t *Qb; ///< Qb[n] S波Q值 real_t *Qainv; ///< 1/Q_p real_t *Qbinv; ///< 1/Q_s + bool *isLiquid; ///< 每层是否为液体层 cplx_t *mu; ///< mu[n] \f$ V_b^2 * \rho \f$ cplx_t *lambda; ///< lambda[n] \f$ V_a^2 * \rho - 2*\mu \f$ diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index 5435052..1acc3c9 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -31,6 +31,7 @@ X(Qb, real_t)\ X(Qainv, real_t)\ X(Qbinv, real_t)\ + X(isLiquid, bool)\ X(mu, cplx_t)\ X(lambda, cplx_t)\ X(delta, cplx_t)\ @@ -209,7 +210,7 @@ void grt_mod1d_xa_xb(GRT_MODEL1D *mod1d, const real_t k) mod1d->caca[i] = caca; mod1d->xa[i] = sqrt(1.0 - caca); - cbcb = (vb > 0.0)? mod1d->c_phase / (vb*atnb) : 0.0; // 考虑液体层 + cbcb = (mod1d->isLiquid[i])? 0.0 : mod1d->c_phase / (vb*atnb); // 考虑液体层 cbcb *= cbcb; mod1d->cbcb[i] = cbcb; mod1d->xb[i] = sqrt(1.0 - cbcb); @@ -390,6 +391,11 @@ GRT_MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, rea mod1d->depsrc = depsrc; mod1d->deprcv = deprcv; + // 记录每层是否为液体层 + for(size_t i = 0; i < mod1d->n; ++i){ + mod1d->isLiquid[i] = (mod1d->Vb[i] == 0.0); + } + // 检查，接收点不能位于液-液、固-液界面 if(ircv < nlay-1 && mod1d->Thk[ircv] == 0.0 && mod1d->Vb[ircv]*mod1d->Vb[ircv+1] == 0.0){ GRTRaiseError( diff --git a/pygrt/C_extension/src/dynamic/grt_greenfn.c b/pygrt/C_extension/src/dynamic/grt_greenfn.c index e8c8dcb..e3f8073 100644 --- a/pygrt/C_extension/src/dynamic/grt_greenfn.c +++ b/pygrt/C_extension/src/dynamic/grt_greenfn.c @@ -985,7 +985,7 @@ int greenfn_main(int argc, char **argv) { // 当震源位于液体层中时，仅允许计算爆炸源对应的格林函数 // 程序结束前会输出对应警告 - if(mod1d->Vb[mod1d->isrc]==0.0){ + if(mod1d->isLiquid[mod1d->isrc]){ Ctrl->G.doHF = Ctrl->G.doVF = Ctrl->G.doDC = false; } @@ -1223,7 +1223,7 @@ int greenfn_main(int argc, char **argv) { } // 输出警告：当震源位于液体层中时，仅允许计算爆炸源对应的格林函数 - if(mod1d->Vb[mod1d->isrc]==0.0){ + if(mod1d->isLiquid[mod1d->isrc]){ GRTRaiseWarning( "The source is located in the liquid layer, " "therefore only the Green's Funtions for the Explosion source will be computed.\n"); diff --git a/pygrt/C_extension/src/dynamic/layer.c b/pygrt/C_extension/src/dynamic/layer.c index 7e31703..9e2bfde 100644 --- a/pygrt/C_extension/src/dynamic/layer.c +++ b/pygrt/C_extension/src/dynamic/layer.c @@ -29,9 +29,9 @@ T N##2 = mod1d->N[iy];\ -static int __freeBound_R(cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML) +static int __freeBound_R(bool isLiquid, cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML) { - if(cbcb0 != 0.0){ + if( ! isLiquid){ // 固体表面 // 公式(5.3.10-14) cplx_t Delta = 0.0; @@ -58,9 +58,9 @@ static int __freeBound_R(cplx_t cbcb0, cplx_t xa0, cplx_t xb0, real_t adsgn, cpl 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) +static int __rigidBound_R(bool isLiquid, cplx_t xa0, cplx_t xb0, real_t adsgn, cplx_t M[2][2], cplx_t *ML) { - if(cbcb0 != 0.0){ + if( ! isLiquid){ // 固体表面 cplx_t Delta = 0.0; @@ -89,11 +89,12 @@ 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]; + bool isLiquid = mod1d->isLiquid[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); + mod1d->M_top.stats = __freeBound_R(isLiquid, 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); + mod1d->M_top.stats = __rigidBound_R(isLiquid, 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); @@ -112,8 +113,9 @@ void grt_botbound_RD(GRT_MODEL1D *mod1d) cplx_t cbcb = mod1d->cbcb[nlay-2]; cplx_t xa = mod1d->xa[nlay-2]; cplx_t xb = mod1d->xb[nlay-2]; + cplx_t isLiquid = mod1d->isLiquid[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); + mod1d->M_bot.stats = __freeBound_R(isLiquid, 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); @@ -149,10 +151,11 @@ void grt_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d) size_t ircv = mod1d->ircv; cplx_t xa = mod1d->xa[ircv]; cplx_t xb = mod1d->xb[ircv]; + bool isLiquid = mod1d->isLiquid[ircv]; real_t k = mod1d->k; cplx_t D11[2][2], D12[2][2]; - if(xb != 1.0){ + if( ! isLiquid){ // 位于固体层 // 公式(5.2.19) D11[0][0] = k; D11[0][1] = k*xb; @@ -181,10 +184,10 @@ void grt_wave2qwv_REV_PSV(GRT_MODEL1D *mod1d) void grt_wave2qwv_REV_SH(GRT_MODEL1D *mod1d) { size_t ircv = mod1d->ircv; - cplx_t xb = mod1d->xb[ircv]; + bool isLiquid = mod1d->isLiquid[ircv]; real_t k = mod1d->k; - if(xb != 1.0){ + if( ! isLiquid){ // 位于固体层 // 公式(5.2.19) if(mod1d->ircvup){// 震源更深 @@ -204,6 +207,7 @@ void grt_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d) size_t ircv = mod1d->ircv; cplx_t xa = mod1d->xa[ircv]; cplx_t xb = mod1d->xb[ircv]; + bool isLiquid = mod1d->isLiquid[ircv]; real_t k = mod1d->k; @@ -217,7 +221,7 @@ void grt_wave2qwv_z_REV_PSV(GRT_MODEL1D *mod1d) cplx_t D12[2][2] = {{-ak, bb}, {aa, -bk}}; // 位于液体层 - if(xb == 1.0){ + if(isLiquid){ D11[0][1] = D11[1][1] = D12[0][1] = D12[1][1] = 0.0; } @@ -236,13 +240,14 @@ void grt_wave2qwv_z_REV_SH(GRT_MODEL1D *mod1d) { size_t ircv = mod1d->ircv; cplx_t xb = mod1d->xb[ircv]; + bool isLiquid = mod1d->isLiquid[ircv]; real_t k = mod1d->k; // 将垂直波函数转为ui,z在(B_m, P_m, C_m)系下的分量 // 新推导的公式 cplx_t bk = k*k*xb; - if(xb != 1.0){ + if( ! isLiquid){ // 位于固体层 if(mod1d->ircvup){// 震源更深 mod1d->uiz_R_EVL = (1.0 - mod1d->M_FA.RUL)*bk; @@ -298,21 +303,21 @@ void grt_RT_matrix_ls_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX * MODEL_2LAYS_ATTRIB(cplx_t, mu); MODEL_2LAYS_ATTRIB(cplx_t, cbcb); MODEL_2LAYS_ATTRIB(real_t, Rho); + MODEL_2LAYS_ATTRIB(bool, isLiquid); // 后缀1表示上层的液体的物理参数，后缀2表示下层的固体的物理参数 // 若mu2==0, 则下层为液体，参数需相互交换 // 讨论液-固 or 固-液 - bool isfluidUp = (mu1 == 0.0); // 上层是否为液体 int sgn = 1; - if(isfluidUp && mu2 == 0.0){ + if(isLiquid1 && isLiquid2){ GRTRaiseError("fluid-fluid interface is not allowed.\n"); } // 使用指针 cplx_t (*pRD)[2], (*pRU)[2]; cplx_t (*pTD)[2], (*pTU)[2]; - if(isfluidUp){ + if(isLiquid1){ pRD = M->RD; pRU = M->RU; pTD = M->TD; pTU = M->TU; } else { @@ -362,23 +367,20 @@ void grt_RT_matrix_ls_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX * void grt_RT_matrix_ls_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { - // TEMPORARY!!!!!! - // 之后不再使用mu或其它变量来判断是否为液体，直接定义一个新的数组 - MODEL_2LAYS_ATTRIB(cplx_t, mu); + MODEL_2LAYS_ATTRIB(bool, isLiquid); // 后缀1表示上层的液体的物理参数，后缀2表示下层的固体的物理参数 // 若mu2==0, 则下层为液体，参数需相互交换 // 讨论液-固 or 固-液 - bool isfluidUp = (mu1 == 0.0); // 上层是否为液体 - if(isfluidUp && mu2 == 0.0){ + if(isLiquid1 && isLiquid2){ GRTRaiseError("fluid-fluid interface is not allowed.\n"); } // 使用指针 cplx_t *pRDL, *pRUL; cplx_t *pTDL, *pTUL; - if(isfluidUp){ + if(isLiquid1){ pRDL = &M->RDL; pRUL = &M->RUL; pTDL = &M->TDL; pTUL = &M->TUL; } else { @@ -494,15 +496,13 @@ void grt_RT_matrix_ss_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M void grt_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { - // TEMPORARY!!!!!! - // 之后不再使用mu或其它变量来判断是否为液体，直接定义一个新的数组 - MODEL_2LAYS_ATTRIB(cplx_t, mu); + MODEL_2LAYS_ATTRIB(bool, isLiquid); // 根据界面两侧的具体情况选择函数 - if(mu1 != 0.0 && mu2 != 0.0){ + if(!isLiquid1 && !isLiquid2){ grt_RT_matrix_ss_PSV(mod1d, iy, M); } - else if(mu1 == 0.0 && mu2 == 0.0){ + else if(isLiquid1 && isLiquid2){ grt_RT_matrix_ll_PSV(mod1d, iy, M); } else{ @@ -513,15 +513,13 @@ void grt_RT_matrix_PSV(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) void grt_RT_matrix_SH(const GRT_MODEL1D *mod1d, const size_t iy, RT_MATRIX *M) { - // TEMPORARY!!!!!! - // 之后不再使用mu或其它变量来判断是否为液体，直接定义一个新的数组 - MODEL_2LAYS_ATTRIB(cplx_t, mu); + MODEL_2LAYS_ATTRIB(bool, isLiquid); // 根据界面两侧的具体情况选择函数 - if(mu1 != 0.0 && mu2 != 0.0){ + if(!isLiquid1 && !isLiquid2){ grt_RT_matrix_ss_SH(mod1d, iy, M); } - else if(mu1 == 0.0 && mu2 == 0.0){ + else if(isLiquid1 && isLiquid2){ grt_RT_matrix_ll_SH(M); } else{ diff --git a/pygrt/c_structures.py b/pygrt/c_structures.py index 8d17400..7c6797d 100755 --- a/pygrt/c_structures.py +++ b/pygrt/c_structures.py @@ -111,6 +111,7 @@ class c_GRT_MODEL1D(Structure): ('Qb', PREAL), ('Qainv', PREAL), ('Qbinv', PREAL), + ('isLiquid', c_bool), ('mu', PCPLX), ('lambda', PCPLX),