diff --git a/docs/source/Gallery/ex15/plot.py b/docs/source/Gallery/ex15/plot.py index aa9609bf..959aff86 100644 --- a/docs/source/Gallery/ex15/plot.py +++ b/docs/source/Gallery/ex15/plot.py @@ -50,8 +50,8 @@ ax.plot(distarr, dynamic_Z, 'k', label='Dynamic Z') ax.plot(distarr, dynamic_R, 'k', label='Dynamic R') -ax.plot(distarr, static_grn['variables']['EXZ']['data'][0] * coef, 'ro', ms=4, label='Static Z') -ax.plot(distarr, static_grn['variables']['EXR']['data'][0] * coef, 'bo', ms=4, label='Static R') +ax.plot(distarr, static_grn['variables']['EXZ']['data'][0, 0, 0] * coef, 'ro', ms=4, label='Static Z') +ax.plot(distarr, static_grn['variables']['EXR']['data'][0, 0, 0] * coef, 'bo', ms=4, label='Static R') ax.set_xlim(0, 50) ax.set_xlabel('Distance (km)') diff --git a/docs/source/Gallery/ex15/plot_all.py b/docs/source/Gallery/ex15/plot_all.py index 855967d9..8f4a2838 100644 --- a/docs/source/Gallery/ex15/plot_all.py +++ b/docs/source/Gallery/ex15/plot_all.py @@ -62,12 +62,12 @@ ms = 2 lw = 0.8 ax.plot(distarr, dynamic_Z, 'k', lw=lw, label='Dynamic Z') - ax.plot(distarr, static_grn['variables'][f'{src}Z']['data'][0] * coef, 'ro', ms=ms, label='Static Z') + ax.plot(distarr, static_grn['variables'][f'{src}Z']['data'][0, 0, 0] * coef, 'ro', ms=ms, label='Static Z') ax.plot(distarr, dynamic_R, 'k', lw=lw, label='Dynamic R') - ax.plot(distarr, static_grn['variables'][f'{src}R']['data'][0] * coef, 'bo', ms=ms, label='Static R') + ax.plot(distarr, static_grn['variables'][f'{src}R']['data'][0, 0, 0] * coef, 'bo', ms=ms, label='Static R') if src not in ['EX', 'VF', 'DD']: ax.plot(distarr, dynamic_T, 'k', lw=lw, label='Dynamic T') - ax.plot(distarr, static_grn['variables'][f'{src}T']['data'][0] * coef, 'go', ms=ms, label='Static T') + ax.plot(distarr, static_grn['variables'][f'{src}T']['data'][0, 0, 0] * coef, 'go', ms=ms, label='Static T') ax.set_xlim(0, 50) ax.set_xlabel('Distance (km)') diff --git a/docs/source/Gallery/ex16/plot.py b/docs/source/Gallery/ex16/plot.py index 84c26fa1..86e0041d 100644 --- a/docs/source/Gallery/ex16/plot.py +++ b/docs/source/Gallery/ex16/plot.py @@ -115,8 +115,8 @@ ax.set_ymargin(0.3) ax = axs2[i] - ax.plot(easts, static1['variables'][chLst[i]]['data'][0], **prop1) - ax.plot(easts, static2['variables'][chLst[i]]['data'][0] * sgn, **prop2) + ax.plot(easts, static1['variables'][chLst[i]]['data'][0, 0, 0], **prop1) + ax.plot(easts, static2['variables'][chLst[i]]['data'][0, 0, 0] * sgn, **prop2) ax.text(0.96, 0.9, chLst[i], transform=ax.transAxes, ha='right', va='top', bbox=dict(fc='w')) ax.ticklabel_format(axis='y', style='sci', scilimits=(0,0)) diff --git a/docs/source/Tutorial/static/run/run.py b/docs/source/Tutorial/static/run/run.py index 30e6ed17..f40dccaa 100755 --- a/docs/source/Tutorial/static/run/run.py +++ b/docs/source/Tutorial/static/run/run.py @@ -24,7 +24,7 @@ print(static_grn.keys()) # dict_keys(['dimensions', 'variables', 'attributes']) print(list(static_grn["variables"].keys())) -# ['north', 'east', 'EXZ', 'EXR', ...] +# ['model', 'depsrc', 'deprcv', 'north', 'east', 'src_va', ..., 'EXZ', 'EXR', ...] # END GRN # --------------------------------------------------------------------------------- diff --git a/docs/source/Tutorial/static/run/run.sh b/docs/source/Tutorial/static/run/run.sh index 4dff5230..24b3f701 100755 --- a/docs/source/Tutorial/static/run/run.sh +++ b/docs/source/Tutorial/static/run/run.sh @@ -118,7 +118,7 @@ grt static syn -S1e24 -M33/50 -N -Gstgrn.nc -Ostsyn_ts.nc gmt begin syn_ts pdf gmtplot_static stsyn_ts.nc -Si0.03c gmt meca -Sz0.5c < top + * 调用方负责 grt_finite_fault_free + * + * @param[in] path 文件路径 + * @param[out] nfault 读入的断层段数 + * @return 新分配的 FINITE_FAULT 数组，失败则报错退出 + */ +FINITE_FAULT *grt_finite_fault_load_coulomb(const char *path, size_t *nfault); + +/** + * 释放 grt_finite_fault_load_coulomb 返回的数组 + * + * @param[in,out] faults 断层数组，可为 NULL + */ +void grt_finite_fault_free(FINITE_FAULT *faults); + +/** + * 由断层几何与 dL/dW 得到倾向/走向长度及剖分数 + * + * 要求 dip ∈ (0, 90]、bot > top，且 dL/dW > 0 + * + * @param[in] fault 有限断层（需已有 dip / 端点坐标） + * @param[in] dL 沿走向剖分间隔 (km) + * @param[in] dW 沿倾向剖分间隔 (km) + * @param[out] W 沿倾向长度 (km) + * @param[out] L 沿走向长度 (km) + * @param[out] nW 倾向方向子断层数 + * @param[out] nL 走向方向子断层数 + */ +void grt_finite_fault_subdiv( + const FINITE_FAULT *fault, real_t dL, real_t dW, + real_t *W, real_t *L, size_t *nW, size_t *nL); + +/** + * 计算 (iW, iL) 子断层中心几何与矩势 + * + * 要求 fault 已 set_derived，且 W/L/nW/nL 由 grt_finite_fault_subdiv 给出 + * 末块可短于 dL/dW，中心取该块中点：i*d + size/2 + * + * @param[in] fault 有限断层 + * @param[in] dL 沿走向剖分间隔 (km) + * @param[in] dW 沿倾向剖分间隔 (km) + * @param[in] W 沿倾向总长 (km) + * @param[in] L 沿走向总长 (km) + * @param[in] iW 倾向方向子断层索引 [0, nW) + * @param[in] iL 走向方向子断层索引 [0, nL) + * @param[out] sub 子断层几何 + */ +void grt_finite_fault_subfault( + const FINITE_FAULT *fault, + real_t dL, real_t dW, real_t W, real_t L, + size_t iW, size_t iL, + FINITE_SUBFAULT *sub); diff --git a/pygrt/C_extension/include/grt/common/model.h b/pygrt/C_extension/include/grt/common/model.h index 2e58ff0a..364e2458 100755 --- a/pygrt/C_extension/include/grt/common/model.h +++ b/pygrt/C_extension/include/grt/common/model.h @@ -128,6 +128,32 @@ void grt_free_mod1d(MODEL1D *mod1d); */ MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid); +/** + * 从模型文件读取 ``nlayer × 6`` 矩阵（Thk, Va, Vb, Rho, Qa, Qb） + * 若首列以层顶深度给出则先转为厚度；末层厚度保持文件原值 + * + * @param[in] modelpath 模型文件路径 + * @param[out] nlayer 层数 + * @param[in] allowLiquid 是否允许液体层（Vs==0） + * @return 新分配的矩阵，调用方负责 free + */ +real_t (* grt_read_modarr_from_file( + const char *modelpath, size_t *nlayer, bool allowLiquid))[GRT_MODARR_NCOL]; + +/** + * 由 ``nlayer × 6`` 模型矩阵按深度查层介质 + * 末层视为半空间；恰落在层界面时取上层 + * + * @param[in] nlayer 层数 + * @param[in] modarr 模型矩阵，每行 Thk/Va/Vb/Rho/Qa/Qb + * @param[in] depth 深度 (km)，须 >= 0 + * @param[out] va P 波速 (km/s)，可为 NULL + * @param[out] vb S 波速 (km/s)，可为 NULL + * @param[out] rho 密度 (g/cm^3)，可为 NULL + */ +void grt_modarr_medium_at_depth( + size_t nlayer, const real_t (*modarr)[GRT_MODARR_NCOL], + real_t depth, real_t *va, real_t *vb, real_t *rho); /** * 设置模型的边界条件，并对底界面做检查 diff --git a/pygrt/C_extension/include/grt/common/search.h b/pygrt/C_extension/include/grt/common/search.h index 2a0e986c..47bfd633 100644 --- a/pygrt/C_extension/include/grt/common/search.h +++ b/pygrt/C_extension/include/grt/common/search.h @@ -157,3 +157,21 @@ typedef int (*grt_compare_fn)(const void *a, const void *b); int grt_argsort( const void *base, size_t n, size_t element_size, grt_compare_fn compare, size_t *indices); + +/** + * 在升序数组上定位 1D 线性插值的左右邻点与权重 + * + * n==1 时要求 |q-x[0]| 在容差内，此时 i0=i1=0、w=0； + * q 越界（超出端点容差）返回 false + * + * @param[in] x 升序坐标数组 + * @param[in] n 长度 + * @param[in] q 查询值 + * @param[out] i0 左邻索引 + * @param[out] i1 右邻索引 + * @param[out] w 权重，使值 = (1-w)*x[i0] + w*x[i1] + * @return 是否在范围内 + */ +bool grt_locateLinearInterp( + const real_t *x, size_t n, real_t q, + size_t *i0, size_t *i1, real_t *w); diff --git a/pygrt/C_extension/include/grt/static/recv_points.h b/pygrt/C_extension/include/grt/static/recv_points.h new file mode 100644 index 00000000..ed7d1d29 --- /dev/null +++ b/pygrt/C_extension/include/grt/static/recv_points.h @@ -0,0 +1,79 @@ +/** + * @file recv_points.h + * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) + * @date 2026-08 + * + * 静态 syn / 后处理用的接收点列表 + * 网格 (-X/-Y 或延用库坐标) 与任意点文件 (-Q) 均展开为点列 + * + */ + +#pragma once + +#include +#include + +#include "grt/common/const.h" + +/** layout 字符串：写入 nc 全局属性，读端据此分支 */ +#define GRT_RECV_LAYOUT_GRID "grid" +#define GRT_RECV_LAYOUT_POINTS "points" + +/** + * 接收点坐标列表（三个分量各自连续存放） + * + * norths/easts/depths 长度均为 npts，单位 km + * is_grid 为真时可由 nnorth/neast 还原二维网格，ipt = ieast + inorth*neast + */ +typedef struct { + size_t npts; + real_t *norths; + real_t *easts; + real_t *depths; + + bool is_grid; + size_t nnorth; + size_t neast; +} GRT_RECV_POINTS; + +/** + * 由 north/east 轴与单一深度展开为点列（is_grid=true） + * + * @param[in] nnorth north 方向点数 + * @param[in] norths north 坐标 (km) + * @param[in] neast east 方向点数 + * @param[in] easts east 坐标 (km) + * @param[in] depth 接收深度 (km) + * @return 新分配的 GRT_RECV_POINTS*，调用方负责 grt_recv_points_free + */ +GRT_RECV_POINTS *grt_recv_points_from_grid( + size_t nnorth, const real_t *norths, + size_t neast, const real_t *easts, + real_t depth); + +/** + * 从 ASCII 文件读任意接收点（is_grid=false） + * + * 每行 north east depth (km)，# 开头为注释 + * + * @param[in] path 文件路径 + * @return 新分配的 GRT_RECV_POINTS* + */ +GRT_RECV_POINTS *grt_recv_points_from_file(const char *path); + +/** + * 释放 GRT_RECV_POINTS（含 norths/easts/depths） + * + * @param[in,out] pts 可为 NULL + */ +void grt_recv_points_free(GRT_RECV_POINTS *pts); + +/** + * 从已打开的 nc 判断是否为 points 布局 + * + * 优先读全局属性 layout；若无属性则看是否存在 point 维 + * + * @param[in] ncid 已打开的 nc id + * @return true 表示 points，false 表示 grid + */ +bool grt_recv_nc_is_points(int ncid); diff --git a/pygrt/C_extension/include/grt/static/static_grn.h b/pygrt/C_extension/include/grt/static/static_grn.h index 986c545a..e3e04bd0 100644 --- a/pygrt/C_extension/include/grt/static/static_grn.h +++ b/pygrt/C_extension/include/grt/static/static_grn.h @@ -41,4 +41,4 @@ void grt_integ_static_grn( realChnlGrid grn[nr], realChnlGrid grn_uiz[nr], realChnlGrid grn_uir[nr], - const char *statsstr); \ No newline at end of file + const char *statsstr); diff --git a/pygrt/C_extension/include/grt/static/stgrnlib.h b/pygrt/C_extension/include/grt/static/stgrnlib.h new file mode 100644 index 00000000..1c5113c3 --- /dev/null +++ b/pygrt/C_extension/include/grt/static/stgrnlib.h @@ -0,0 +1,133 @@ +/** + * @file stgrnlib.h + * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) + * @date 2026-08 + * + * 静态格林函数库 STGRNLIB + * + * nc 维度为 depsrc×deprcv×north×east； + * C 侧深度/水平坐标数组为 depsrcs、deprcvs、norths、easts； + * 物理上格林函数还依赖震中距 r = hypot(north, east)； + * 用 -R 建库时通常为 nnorth=1, norths=[0], easts=R（震中距序列） + * + */ + +#pragma once + +#include +#include "grt/common/const.h" + +/** + * 静态格林函数库 + * + * 数组布局： + * - depsrcs[ndepsrc], deprcvs[ndeprcv] 严格升序 (km) + * - norths[nnorth], easts[neast] 为 north / east 坐标 (km)，各自严格升序 + * - rs[nr] 为网格序震中距，ipt = ieast + inorth*neast，nr = nnorth*neast + * - sort_rs / sort_rs_idx 为 rs 的升序排列及回指网格 ipt 的索引 + * - isUniform / dr：网格序 rs 是否已是等距升序及其步长（同 syn 查找加速） + * - u[is][ir][ipt][im][c] + * - 分量符号约定与 nc 文件 / Python dict 一致（Z 已取反等） + */ +typedef struct { + size_t ndepsrc; ///< 震源深度数量 + real_t *depsrcs; ///< depsrcs[ndepsrc] (km) + + size_t ndeprcv; ///< 接收深度数量 + real_t *deprcvs; ///< deprcvs[ndeprcv] (km) + + size_t nnorth; ///< north 方向点数 + size_t neast; ///< east 方向点数（-R 建库时即震中距点数） + real_t *norths; ///< north 坐标 norths[nnorth] + real_t *easts; ///< east 坐标 easts[neast]（-R 建库时即震中距） + + size_t nr; ///< = nnorth * neast + real_t *rs; ///< 网格序震中距 rs[nr] + real_t *sort_rs; ///< 升序震中距 sort_rs[nr] + size_t *sort_rs_idx; ///< sort_rs[i] 对应网格 ipt = sort_rs_idx[i] + bool isUniform; ///< 网格序 rs 是否已是等距升序 + real_t dr; ///< 等距步长；非等距时无意义 + + bool calc_upar; ///< 是否含位移空间导数 + + real_t *src_va; ///< 各震源深度处 P 波速 src_va[ndepsrc] + real_t *src_vb; ///< 各震源深度处 S 波速 + real_t *src_rho; ///< 各震源深度处密度 + + real_t *rcv_va; ///< 各接收深度处 P 波速 rcv_va[ndeprcv] + real_t *rcv_vb; ///< 各接收深度处 S 波速 + real_t *rcv_rho; ///< 各接收深度处密度 + + size_t nlayer; ///< 建库模型层数 + real_t (*modarr)[GRT_MODARR_NCOL]; ///< 模型矩阵 [nlayer][6]：Thk/Va/Vb/Rho/Qa/Qb + + realChnlGrid ***u; ///< u[is][ir]，每个为 realChnlGrid[nr] + realChnlGrid ***uiz; ///< 可选，calc_upar=false 时为 NULL + realChnlGrid ***uir; ///< 可选，calc_upar=false 时为 NULL +} STGRNLIB; + + +/** + * 按维度申请空壳 STGRNLIB：拷贝坐标轴，介质数组清零，并分配 u/uiz/uir + * + * 调用方随后填入介质与格林函数；最终用 grt_stgrnlib_free 释放 + * + * @param[in] ndepsrc 震源深度点数 + * @param[in] depsrcs 震源深度数组 (km) + * @param[in] ndeprcv 接收深度点数 + * @param[in] deprcvs 接收深度数组 (km) + * @param[in] nnorth north 方向点数 + * @param[in] norths north 坐标数组 (km) + * @param[in] neast east 方向点数 + * @param[in] easts east 坐标数组 (km) + * @param[in] calc_upar 是否分配位移偏导 + * @return 新分配的 STGRNLIB* + */ +STGRNLIB *grt_stgrnlib_alloc( + size_t ndepsrc, const real_t *depsrcs, + size_t ndeprcv, const real_t *deprcvs, + size_t nnorth, const real_t *norths, + size_t neast, const real_t *easts, + bool calc_upar); + +/** + * 释放 STGRNLIB 内部所有堆内存，并 free(lib) 本身 + * （适用于 load_nc / alloc 返回的指针） + * + * @param[in,out] lib STGRNLIB 结构体，可为 NULL + */ +void grt_stgrnlib_free(STGRNLIB *lib); + +/** + * 设置建库模型矩阵（会拷贝一份） + * + * @param[in,out] lib STGRNLIB + * @param[in] nlayer 层数，须 > 0 + * @param[in] modarr 模型矩阵，每行 Thk/Va/Vb/Rho/Qa/Qb + */ +void grt_stgrnlib_set_modarr( + STGRNLIB *lib, size_t nlayer, const real_t (*modarr)[GRT_MODARR_NCOL]); + +/** + * 由震中距与震源深度采样推断默认子断层尺寸 min(dr, dz) + * + * @param[in] lib STGRNLIB + * @return 默认 dL=dW (km) + */ +real_t grt_stgrnlib_default_subfault_size(const STGRNLIB *lib); + +/** + * 从四维 nc 文件加载完整 STGRNLIB + * + * @param[in] path nc 文件路径 + * @return 新分配的 STGRNLIB*，调用方负责 grt_stgrnlib_free + */ +STGRNLIB *grt_stgrnlib_load_nc(const char *path); + +/** + * 将 STGRNLIB 写为四维 nc 文件 + * + * @param[in] lib 已填充的库 + * @param[in] path 输出路径 + */ +void grt_stgrnlib_save_nc(const STGRNLIB *lib, const char *path); diff --git a/pygrt/C_extension/src/common/finite_fault.c b/pygrt/C_extension/src/common/finite_fault.c new file mode 100644 index 00000000..08b80526 --- /dev/null +++ b/pygrt/C_extension/src/common/finite_fault.c @@ -0,0 +1,152 @@ +/** + * @file finite_fault.c + * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) + * @date 2026-08 + * + * Coulomb 格式有限断层：读入、衍生量与几何剖分 + * + */ + +#include +#include +#include +#include + +#include "grt/common/finite_fault.h" +#include "grt/common/util.h" + +/** 检查 dip ∈ (0, 90] 且 bot > top */ +static void check_fault_geometry(const FINITE_FAULT *f, const char *where) +{ + if(f->dip <= 0.0 || f->dip > 90.0){ + GRTRaiseError("%s: dip (%.6g deg) must be in (0, 90].", where, f->dip); + } + if(!(f->bot > f->top)){ + GRTRaiseError("%s: bot (%.6g km) must be greater than top (%.6g km).", where, f->bot, f->top); + } +} + + +void grt_finite_fault_set_derived(FINITE_FAULT *f) +{ + f->strike = 1.0 / DEG1 * atan2(f->east_end - f->east_begin, f->north_end - f->north_begin); + f->rake = 1.0 / DEG1 * atan2(f->reverse, - f->right_lateral); + f->slip = hypot(f->right_lateral, f->reverse); // m +} + + +FINITE_FAULT *grt_finite_fault_load_coulomb(const char *path, size_t *nfault) +{ + if(path == NULL || nfault == NULL){ + GRTRaiseError("path/nfault is NULL."); + } + + GRTCheckFileExist(path); + FILE *fp = GRTCheckOpenFile(path, "r"); + + char *line = NULL; + size_t nlen = 0; + // 跳过两行表头（列名 + 占位行） + for(int ih = 0; ih < 2; ++ih){ + if(grt_getline(&line, &nlen, fp) <= 0){ + fclose(fp); + GRT_SAFE_FREE_PTR(line); + GRTRaiseError("read header of %s failed.", path); + } + } + + FINITE_FAULT *faults = NULL; + size_t n = 0; + while(grt_getline(&line, &nlen, fp) != -1){ + faults = (FINITE_FAULT *)realloc(faults, sizeof(FINITE_FAULT) * (n + 1)); + FINITE_FAULT *f = faults + n; + memset(f, 0, sizeof(*f)); + + real_t dum1, dum2; + int nscan = sscanf(line, "%lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf", + &dum1, &f->east_begin, &f->north_begin, &f->east_end, &f->north_end, + &dum2, &f->right_lateral, &f->reverse, &f->dip, &f->top, &f->bot); + if(nscan != 11){ + fclose(fp); + GRT_SAFE_FREE_PTR(line); + GRT_SAFE_FREE_PTR(faults); + GRTRaiseError("parse line %zu of %s failed.", n + 3, path); + } + + char where[256]; + snprintf(where, sizeof(where), "in %s line %zu", path, n + 3); + check_fault_geometry(f, where); + + grt_finite_fault_set_derived(f); + n++; + } + + GRT_SAFE_FREE_PTR(line); + fclose(fp); + + if(n == 0){ + GRTRaiseError("no fault in %s.", path); + } + + *nfault = n; + return faults; +} + + +void grt_finite_fault_free(FINITE_FAULT *faults) +{ + free(faults); +} + + +void grt_finite_fault_subdiv( + const FINITE_FAULT *fault, real_t dL, real_t dW, + real_t *W, real_t *L, size_t *nW, size_t *nL) +{ + if(dL <= 0.0 || dW <= 0.0){ + GRTRaiseError("dL and dW must be positive."); + } + check_fault_geometry(fault, "finite fault"); + + *W = (fault->bot - fault->top) / sin(DEG1 * fault->dip); + *L = hypot(fault->east_end - fault->east_begin, fault->north_end - fault->north_begin); + if(*W <= 0.0 || *L <= 0.0){ + GRTRaiseError("fault along-dip/along-strike length must be positive (W=%.6g, L=%.6g).", *W, *L); + } + *nW = GRT_MAX(1, (size_t)ceil(*W / dW)); + *nL = GRT_MAX(1, (size_t)ceil(*L / dL)); +} + + +void grt_finite_fault_subfault( + const FINITE_FAULT *fault, + real_t dL, real_t dW, real_t W, real_t L, + size_t iW, size_t iL, + FINITE_SUBFAULT *sub) +{ + // 末块可短于 dW/dL，中心取该块中点：i*d + size/2，而不是 (i+0.5)*size + real_t width = GRT_MIN(dW, W - iW * dW); + real_t length = GRT_MIN(dL, L - iL * dL); + if(width <= 0.0 || length <= 0.0){ + GRTRaiseError("nonpositive subfault size at (iW=%zu, iL=%zu).", iW, iL); + } + real_t w = iW * dW + 0.5 * width; + real_t l = iL * dL + 0.5 * length; + + real_t sind = sin(DEG1 * fault->dip); + real_t cosd = cos(DEG1 * fault->dip); + real_t sins = sin(DEG1 * fault->strike); + real_t coss = cos(DEG1 * fault->strike); + + real_t hproj = w * cosd; + real_t east0 = fault->east_begin + l * sins; + real_t north0 = fault->north_begin + l * coss; + + sub->width = width; + sub->length = length; + sub->depsrc = fault->top + w * sind; + sub->east = east0 + hproj * coss; + sub->north = north0 - hproj * sins; + // slip(m) * width(km) * length(km) → cm^3 + sub->potency = fault->slip * width * length * 1e12; +} diff --git a/pygrt/C_extension/src/common/model.c b/pygrt/C_extension/src/common/model.c index 781ae291..577838f4 100755 --- a/pygrt/C_extension/src/common/model.c +++ b/pygrt/C_extension/src/common/model.c @@ -86,6 +86,7 @@ void grt_realloc_mod1d(MODEL1D *mod1d, size_t n) void grt_free_mod1d(MODEL1D *mod1d) { + if(mod1d == NULL) return; #define X(P, T) GRT_SAFE_FREE_PTR(mod1d->P); __MODEL1D_FOR_EACH_ARRAY #undef X @@ -94,63 +95,53 @@ void grt_free_mod1d(MODEL1D *mod1d) } -MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t deprcv, bool allowLiquid) +real_t (* grt_read_modarr_from_file( + const char *modelpath, size_t *nlayer, bool allowLiquid))[GRT_MODARR_NCOL] { GRTCheckFileExist(modelpath); - - if(depsrc * deprcv < 0.0){ - GRTRaiseError("depsrc and deprcv should have the same sign."); + if(nlayer == NULL){ + GRTRaiseError("nlayer is NULL."); } - - FILE *fp = GRTCheckOpenFile(modelpath, "r"); - - // 初始化 - MODEL1D *mod1d = grt_init_mod1d(1); + FILE *fp = GRTCheckOpenFile(modelpath, "r"); - const int ncols = 6; // 模型文件有6列，或除去qa qb有四列 - const int ncols_noQ = 4; + const int ncols = GRT_MODARR_NCOL; + const int ncols_noQ = GRT_MODARR_NCOL - 2; // 不含 Qa, Qb size_t iline = 0; - real_t h, va, vb, rho, qa, qb; - real_t (*modarr)[ncols] = NULL; - h = va = vb = rho = qa = qb = 0.0; size_t nlay = 0; - mod1d->io_depth = false; + bool io_depth = false; + real_t (*modarr)[GRT_MODARR_NCOL] = NULL; + real_t h, va, vb, rho, qa, qb; - size_t len; + size_t len = 0; char *line = NULL; - - while(grt_getline(&line, &len, fp) != -1) { + while(grt_getline(&line, &len, fp) != -1){ iline++; - - // 注释行 - if(grt_is_comment_or_empty(line)) continue; + if(grt_is_comment_or_empty(line)) continue; h = va = vb = rho = qa = qb = 0.0; int nscan = sscanf(line, "%lf %lf %lf %lf %lf %lf\n", &h, &va, &vb, &rho, &qa, &qb); if(ncols != nscan && ncols_noQ != nscan){ GRTRaiseError("Model file read error in line %zu.\n", iline); - }; + } - // 读取首行，如果首行首列为 0 ，则首列指示每层顶界面深度而非厚度 + // 首行首列为 0 时，首列表示层顶深度而非厚度 if(nlay == 0 && h == 0.0){ - mod1d->io_depth = true; + io_depth = true; } if(va <= 0.0 || rho <= 0.0 || (ncols == nscan && (qa <= 0.0 || qb <= 0.0))){ GRTRaiseError("In model file, line %zu, nonpositive value is not supported.\n", iline); } - if(vb < 0.0){ GRTRaiseError("In model file, line %zu, negative Vs is not supported.\n", iline); } - if(!allowLiquid && vb == 0.0){ GRTRaiseError("In model file, line %zu, Vs==0.0 is not supported.\n", iline); } - modarr = (real_t(*)[ncols])realloc(modarr, sizeof(real_t)*ncols*(nlay+1)); - + modarr = (real_t (*)[GRT_MODARR_NCOL])realloc( + modarr, sizeof(real_t) * GRT_MODARR_NCOL * (nlay + 1)); modarr[nlay][0] = h; modarr[nlay][1] = va; modarr[nlay][2] = vb; @@ -158,26 +149,65 @@ MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t modarr[nlay][4] = qa; modarr[nlay][5] = qb; nlay++; - } + fclose(fp); + GRT_SAFE_FREE_PTR(line); - if(iline==0 || modarr==NULL){ + if(iline == 0 || modarr == NULL){ GRTRaiseError("Model file %s read error.\n", modelpath); } - // 如果读取了深度，转为厚度 - if(mod1d->io_depth){ - for(size_t i=1; iio_depth = false; // 已在 read_modarr 中转为厚度 + + real_t h, va, vb, rho, qa, qb; size_t isrc=0, ircv=0; size_t *pmin_idx, *pmax_idx, *pimg_idx; real_t depth = 0.0, depmin, depmax, depimg; @@ -293,9 +323,7 @@ MODEL1D * grt_read_mod1d_from_file(const char *modelpath, real_t depsrc, real_t depth += mod1d->Thk[iz]; } - fclose(fp); GRT_SAFE_FREE_PTR(modarr); - GRT_SAFE_FREE_PTR(line); // 设置一个默认边界条件 mod1d->topbound = GRT_BOUND_FREE; diff --git a/pygrt/C_extension/src/common/search.c b/pygrt/C_extension/src/common/search.c index 5d939850..fabbdc32 100644 --- a/pygrt/C_extension/src/common/search.c +++ b/pygrt/C_extension/src/common/search.c @@ -9,6 +9,7 @@ #include #include #include +#include #include "grt/common/search.h" #include "grt/common/const.h" @@ -227,3 +228,43 @@ int grt_argsort( free(pairs); return 0; } + + +bool grt_locateLinearInterp( + const real_t *x, size_t n, real_t q, + size_t *i0, size_t *i1, real_t *w) +{ + const real_t atol = 1e-8; + + if(x == NULL || n == 0 || i0 == NULL || i1 == NULL || w == NULL){ + return false; + } + if(n == 1){ + if(fabs(q - x[0]) > atol) return false; + *i0 = *i1 = 0; + *w = 0.0; + return true; + } + if(q < x[0] - atol || q > x[n - 1] + atol){ + return false; + } + if(q <= x[0]){ + *i0 = *i1 = 0; + *w = 0.0; + return true; + } + if(q >= x[n - 1]){ + *i0 = *i1 = n - 1; + *w = 0.0; + return true; + } + size_t i = 0; + for(; i + 1 < n; ++i){ + if(q <= x[i + 1] + atol) break; + } + *i0 = i; + *i1 = i + 1; + real_t dx = x[*i1] - x[*i0]; + *w = (fabs(dx) < atol) ? 0.0 : (q - x[*i0]) / dx; + return true; +} diff --git a/pygrt/C_extension/src/static/grt_static_greenfn.c b/pygrt/C_extension/src/static/grt_static_greenfn.c index a2d5aa26..73c2971d 100644 --- a/pygrt/C_extension/src/static/grt_static_greenfn.c +++ b/pygrt/C_extension/src/static/grt_static_greenfn.c @@ -20,16 +20,16 @@ typedef struct { struct { bool active; char *s_modelpath; ///< 模型路径 - const char *s_modelname; ///< 模型名称 - MODEL1D *mod1d; ///< 模型结构体指针 } M; - /** 震源和接收器深度 */ + /** 震源和接收器深度：-Dsrc/rcv 或 -Ds/-Dr */ struct { - bool active; - real_t depsrc; - real_t deprcv; - char *s_depsrc; - char *s_deprcv; + bool active; ///< 旧式 -D/ + bool s_active; ///< -Ds + bool r_active; ///< -Dr + size_t ndepsrc; + real_t *depsrcs; + size_t ndeprcv; + real_t *deprcvs; } D; /** 顶层和底层的边界条件 */ struct { @@ -99,11 +99,10 @@ typedef struct { static void free_Ctrl(GRT_MODULE_CTRL *Ctrl){ // M GRT_SAFE_FREE_PTR(Ctrl->M.s_modelpath); - grt_free_mod1d(Ctrl->M.mod1d); - + // D - GRT_SAFE_FREE_PTR(Ctrl->D.s_depsrc); - GRT_SAFE_FREE_PTR(Ctrl->D.s_deprcv); + GRT_SAFE_FREE_PTR(Ctrl->D.depsrcs); + GRT_SAFE_FREE_PTR(Ctrl->D.deprcvs); // X GRT_SAFE_FREE_PTR(Ctrl->X.norths); @@ -117,9 +116,7 @@ static void free_Ctrl(GRT_MODULE_CTRL *Ctrl){ GRT_SAFE_FREE_PTR(Ctrl->rs); // S - if(Ctrl->S.active){ - GRT_SAFE_FREE_PTR(Ctrl->S.s_statsdir); - } + GRT_SAFE_FREE_PTR(Ctrl->S.s_statsdir); GRT_SAFE_FREE_PTR(Ctrl); } @@ -138,12 +135,16 @@ printf("\n" "\n\n" "Usage:\n" "----------------------------------------------------------------\n" -" grt static greenfn -M -D/ -O \n" +" grt static greenfn -M -O\n" +" (-D/ | -Ds -Dr) \n" " [-X//] [-Y//] \n" " [-R,[,...]|//|]\n" " [-L] [-C[d|p|n]] [-Bf|F|r|R|h|H] \n" " [-K[+k][+f][+e]] [-S] [-e]\n" "\n" +" Output is always one nc file with dims\n" +" [depsrc][deprcv][north][east] (even if each depth size is 1).\n" +"\n" " There're two ways to define the \"epicentral distances\":\n" " 1. set both -X and -Y (north/east). The kernel still depends\n" " only on r=hypot(north,east); results are stored with north/east\n" @@ -167,8 +168,18 @@ printf("\n" " The number of layers are unlimited.\n" "\n" " -D/\n" -" : source depth (km).\n" -" : receiver depth (km).\n" +" Single source/receiver depth (km). Compatible form.\n" +" Mutually exclusive with -Ds/-Dr.\n" +"\n" +" -Ds Source depth list (km), same syntax as -R:\n" +" + z1,z2[,...]\n" +" + z1/z2/dz\n" +" + \n" +" Must be paired with -Dr.\n" +"\n" +" -Dr\n" +" Receiver depth list (km), same syntax as -Ds.\n" +" Must be paired with -Ds.\n" "\n" " -X//\n" " Set the equidistant points in the north direction.\n" @@ -184,9 +195,10 @@ printf("\n" "\n" " -R,[,...]|//|\n" " Multiple epicentral distances (km), support three ways:\n" -" + ,[,...]: seperated by comma.\n" +" + ,[,...]: seperated by comma (strictly ascending).\n" " + //: equal distance within [r1,r2].\n" -" + : each line contains a distance value.\n" +" + : each line contains a distance value\n" +" (must be strictly ascending).\n" "\n" " -O Filepath to output nc grid.\n" "\n" @@ -236,6 +248,7 @@ printf("\n" " Default 0.0 not use.\n" "\n" " -S Output statsfile in wavenumber integration.\n" +" Only available for a single source/receiver depth.\n" "\n" " -e Compute the spatial derivatives, ui_z and ui_r,\n" " of displacement u. In columns, prefix \"r\" means \n" @@ -251,11 +264,96 @@ printf("\n" "----------------------------------------------------------------\n" " grt static greenfn -Mmilrow -D2/0 -X-10/10/1 -Y-10/10/1 -Ostgrn.nc\n" " grt static greenfn -Mmilrow -D2/0 -R0/20/1 -Ostgrn.nc\n" +" grt static greenfn -Mmilrow -Ds0.2/50/0.5 -Dr0 -R0/500/0.5 -Ostgrn.nc -e\n" "\n\n\n" ); } +/** + * 解析深度列表（语法同 -R），结果升序去重写入 *zs / *nz + * + * 三种输入形式按优先级依次尝试： + * 1. 仅含数字与分隔符时：按逗号拆成离散列表 z1,z2,... + * 2. 可扫成 z1/z2/dz：按等间距生成 [z1, z1+dz, ..., <=z2] + * 3. 否则当作文件路径：逐行读入数值 + * + * 随后转为 real_t、禁止负深度、升序排序并按 1e-8 容差去重 + * + * @param[in] optarg -Ds/-Dr 后的字符串（不含前缀 s/r） + * @param[out] zs 新分配的深度数组，调用方释放 + * @param[out] nz 去重后的点数 + * @param[in] optname 用于报错的选项名（'s' 或 'r'） + */ +static void parse_depth_spec(const char *optarg, real_t **zs, size_t *nz, char optname) +{ + real_t a1, a2, delta; + char **s_vals = NULL; + size_t n = 0; + + // 形式 1：逗号分隔列表（字符集与 -R 一致） + if(grt_string_composed_of(optarg, GRT_NUM_STR "eE+-" ".,")){ + s_vals = grt_string_split(optarg, ",", &n); + } + // 形式 2：等间距 z1/z2/dz + else if(3 == sscanf(optarg, "%lf/%lf/%lf", &a1, &a2, &delta)){ + if(delta <= 0){ + GRTRaiseError("-%c: nonpositive spacing (%f).", optname, delta); + } + if(a1 > a2){ + GRTRaiseError("-%c: start (%f) > end (%f).", optname, a1, a2); + } + n = (size_t)floor((a2 - a1) / delta) + 1; + s_vals = (char **)calloc(n, sizeof(char *)); + for(size_t i = 0; i < n; ++i){ + GRT_SAFE_ASPRINTF(&s_vals[i], "%.*f", 8, a1 + delta * i); + } + } + // 形式 3：从文件逐行读取 + else { + FILE *fp = GRTCheckOpenFile(optarg, "r"); + s_vals = grt_string_from_file(fp, &n); + fclose(fp); + } + + if(n == 0){ + GRTRaiseError("-%c: empty depth list.", optname); + } + + // 字符串 -> 数值，并检查非负 + real_t *raw = (real_t *)calloc(n, sizeof(real_t)); + for(size_t i = 0; i < n; ++i){ + raw[i] = atof(s_vals[i]); + if(raw[i] < 0.0){ + GRTRaiseError("-%c: negative depth (%f) is not supported.", optname, raw[i]); + } + } + GRT_SAFE_FREE_PTR_ARRAY(s_vals, n); + + // 升序排序（argsort 写索引，再按索引取数） + size_t *order = (size_t *)calloc(n, sizeof(size_t)); + for(size_t i = 0; i < n; ++i) order[i] = i; + if(n > 1 && grt_argsort(raw, n, sizeof(*raw), grt_compare_real_t, order) != 0){ + GRTRaiseError("-%c: unable to sort depths.", optname); + } + + // 容差去重，保留升序唯一深度 + real_t *uniq = (real_t *)calloc(n, sizeof(real_t)); + size_t nu = 0; + for(size_t i = 0; i < n; ++i){ + real_t v = raw[order[i]]; + if(nu == 0 || fabs(v - uniq[nu - 1]) > 1e-8){ + uniq[nu++] = v; + } + } + GRT_SAFE_FREE_PTR(raw); + GRT_SAFE_FREE_PTR(order); + + *zs = uniq; + *nz = nu; +} + + @@ -277,25 +375,31 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ case 'M': Ctrl->M.active = true; Ctrl->M.s_modelpath = strdup(optarg); - Ctrl->M.s_modelname = grt_get_basename(Ctrl->M.s_modelpath); break; - // 震源和场点深度， -Ddepsrc/deprcv + // -Dsrc/rcv 或 -Ds / -Dr case 'D': - Ctrl->D.active = true; - Ctrl->D.s_depsrc = (char*)malloc(sizeof(char)*(strlen(optarg)+1)); - Ctrl->D.s_deprcv = (char*)malloc(sizeof(char)*(strlen(optarg)+1)); - if(2 != sscanf(optarg, "%[^/]/%s", Ctrl->D.s_depsrc, Ctrl->D.s_deprcv)){ - GRTBadOptionError(D, ""); - }; - if(1 != sscanf(Ctrl->D.s_depsrc, "%lf", &Ctrl->D.depsrc)){ - GRTBadOptionError(D, ""); - } - if(1 != sscanf(Ctrl->D.s_deprcv, "%lf", &Ctrl->D.deprcv)){ - GRTBadOptionError(D, ""); - } - if(Ctrl->D.depsrc < 0.0 || Ctrl->D.deprcv < 0.0){ - GRTBadOptionError(D, "Negative value in -D is not supported."); + if(optarg[0] == 's'){ + Ctrl->D.s_active = true; + parse_depth_spec(optarg + 1, &Ctrl->D.depsrcs, &Ctrl->D.ndepsrc, 's'); + } else if(optarg[0] == 'r'){ + Ctrl->D.r_active = true; + parse_depth_spec(optarg + 1, &Ctrl->D.deprcvs, &Ctrl->D.ndeprcv, 'r'); + } else { + Ctrl->D.active = true; + real_t depsrc, deprcv; + if(2 != sscanf(optarg, "%lf/%lf", &depsrc, &deprcv)){ + GRTBadOptionError(D, ""); + } + if(depsrc < 0.0 || deprcv < 0.0){ + GRTBadOptionError(D, "Negative value in -D is not supported."); + } + Ctrl->D.ndepsrc = 1; + Ctrl->D.ndeprcv = 1; + Ctrl->D.depsrcs = (real_t *)calloc(1, sizeof(real_t)); + Ctrl->D.deprcvs = (real_t *)calloc(1, sizeof(real_t)); + Ctrl->D.depsrcs[0] = depsrc; + Ctrl->D.deprcvs[0] = deprcv; } break; @@ -520,6 +624,12 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } GRT_SAFE_FREE_PTR_ARRAY(s_easts, Ctrl->Y.neast); + for(size_t i = 1; i < Ctrl->Y.neast; ++i){ + if(!(Ctrl->Y.easts[i] > Ctrl->Y.easts[i - 1])){ + GRTBadOptionError(R, "Epicentral distances must be strictly ascending."); + } + } + Ctrl->X.nnorth = 1; Ctrl->X.norths = (real_t*)calloc(Ctrl->X.nnorth, sizeof(real_t)); Ctrl->X.norths[0] = 0.0; @@ -549,11 +659,19 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ // 检查必须设置的参数是否有设置 GRTCheckOptionSet(argc > 1); GRTCheckOptionActive(Ctrl, M); - GRTCheckOptionActive(Ctrl, D); GRTCheckOptionActive(Ctrl, X); GRTCheckOptionActive(Ctrl, Y); GRTCheckOptionActive(Ctrl, O); + // 深度选项：-D 与 -Ds/-Dr 互斥；-Ds/-Dr 必须成对 + if(Ctrl->D.active && (Ctrl->D.s_active || Ctrl->D.r_active)){ + GRTRaiseError("Options -D and -Ds/-Dr are mutually exclusive."); + } else if(Ctrl->D.s_active != Ctrl->D.r_active){ + GRTRaiseError("Options -Ds and -Dr must be set together."); + } else if(!Ctrl->D.active && !Ctrl->D.s_active){ + GRTRaiseError("Depth option required: -D/ or -Ds... -Dr..."); + } + // 设置震中距数组 Ctrl->nr = Ctrl->X.nnorth*Ctrl->Y.neast; Ctrl->rs = (real_t*)calloc(Ctrl->nr, sizeof(real_t)); @@ -567,15 +685,14 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ /** - * 静态积分前准备：默认 Length，填充 K_INTEG_PROCESS - * 不分配输出缓冲，不处理 stats 路径 + * 静态积分前准备：按震中距与用户参数填充深度无关的 K_INTEG_PROCESS 字段 + * 不写入 Kproc->k0；调用方按 hs 自行缩放 */ static void prepare_static_grn( - MODEL1D *mod1d, size_t nr, real_t *rs, real_t Length, real_t filonLength, real_t safilonTol, real_t filonCut, - real_t k0, real_t keps, bool use_kmax_ref, + real_t keps, bool use_kmax_ref, int convmet, K_INTEG_PROCESS *Kproc) { @@ -594,8 +711,6 @@ static void prepare_static_grn( memset(Kproc, 0, sizeof(*Kproc)); { - real_t hs = GRT_MAX(fabs(mod1d->depsrc - mod1d->deprcv), GRT_MIN_DEPTH_GAP_SRC_RCV); - Kproc->k0 = k0 * PI / hs; Kproc->use_kmax_ref = use_kmax_ref; // 显式收敛方法时不使用 keps Kproc->keps = (convmet != K_INTEG_CONVERG_AUTO) ? 0.0 : keps; @@ -617,167 +732,177 @@ static void prepare_static_grn( } -/** 子模块主函数 */ -int static_greenfn_main(int argc, char **argv){ - GRT_MODULE_CTRL *Ctrl = calloc(1, sizeof(*Ctrl)); +/** 将积分结果按符号约定写入 STGRNLIB 的一层 */ +static void copy_grn_slice_with_sign( + STGRNLIB *lib, size_t is, size_t ir, + size_t nr, bool calc_upar, + const realChnlGrid *grn, const realChnlGrid *grn_uiz, const realChnlGrid *grn_uir) +{ + GRT_LOOP_ChnlGrid(im, c){ + int modr = GRT_SRC_M_ORDERS[im]; + if(modr == 0 && GRT_ZRT_CODES[c] == 'T') continue; + int sgn0 = (GRT_ZRT_CODES[c] == 'Z') ? -1 : 1; + for(size_t ipt = 0; ipt < nr; ++ipt){ + lib->u[is][ir][ipt][im][c] = sgn0 * grn[ipt][im][c]; + if(calc_upar){ + lib->uiz[is][ir][ipt][im][c] = (-1) * sgn0 * grn_uiz[ipt][im][c]; + lib->uir[is][ir][ipt][im][c] = sgn0 * grn_uir[ipt][im][c]; + } + } + } +} - // 传入参数 - getopt_from_command(Ctrl, argc, argv); - // 读入模型文件（暂先不考虑液体层） - if((Ctrl->M.mod1d = grt_read_mod1d_from_file(Ctrl->M.s_modelpath, Ctrl->D.depsrc, Ctrl->D.deprcv, false)) == NULL){ - exit(EXIT_FAILURE); +/** + * 循环计算多震源/接收深度静态格林函数并写入单个四维 nc + * + * Kproc 须已由 prepare_static_grn 填好深度无关字段； + * 循环内按各 (depsrc, deprcv) 的 hs 更新局部拷贝的 k0 + */ +static void compute_stgrnlib_to_nc( + const char *modelpath, + size_t ndepsrc, const real_t *depsrcs, + size_t ndeprcv, const real_t *deprcvs, + size_t nnorth, const real_t *norths, + size_t neast, const real_t *easts, + real_t k0, + K_INTEG_PROCESS *Kproc, + GRT_BOUND_TYPE topbound, GRT_BOUND_TYPE botbound, + bool calc_upar, + const char *outpath, + const char *statsstr) +{ + if(modelpath == NULL || outpath == NULL || Kproc == NULL + || depsrcs == NULL || deprcvs == NULL || norths == NULL || easts == NULL){ + GRTRaiseError("NULL argument."); + } + if(ndepsrc == 0 || ndeprcv == 0 || nnorth == 0 || neast == 0){ + GRTRaiseError("empty dimension."); + } + if(statsstr != NULL && (ndepsrc > 1 || ndeprcv > 1)){ + GRTRaiseError("-S / statsstr is only available for a single source/receiver depth."); } - MODEL1D *mod1d = Ctrl->M.mod1d; - // 边界条件 - grt_set_mod1d_boundary(mod1d, Ctrl->B.topbound, Ctrl->B.botbound); + STGRNLIB *lib = grt_stgrnlib_alloc( + ndepsrc, depsrcs, ndeprcv, deprcvs, nnorth, norths, neast, easts, calc_upar); + size_t nr = lib->nr; + real_t *rs = lib->rs; - // 波数积分输出目录 - if(Ctrl->S.active){ - Ctrl->S.s_statsdir = NULL; - GRT_SAFE_ASPRINTF(&Ctrl->S.s_statsdir, "stgrtstats"); - // 建立保存目录 - GRTCheckMakeDir(Ctrl->S.s_statsdir); - GRT_SAFE_ASPRINTF(&Ctrl->S.s_statsdir, "%s/%s_%s_%s", Ctrl->S.s_statsdir, Ctrl->M.s_modelname, Ctrl->D.s_depsrc, Ctrl->D.s_deprcv); - GRTCheckMakeDir(Ctrl->S.s_statsdir); - } + realChnlGrid *grn = (realChnlGrid *)calloc(nr, sizeof(*grn)); + realChnlGrid *grn_uiz = calc_upar ? (realChnlGrid *)calloc(nr, sizeof(*grn_uiz)) : NULL; + realChnlGrid *grn_uir = calc_upar ? (realChnlGrid *)calloc(nr, sizeof(*grn_uir)) : NULL; - // 建立格林函数的浮点数 - realChnlGrid *grn = (realChnlGrid *) calloc(Ctrl->nr, sizeof(*grn)); - realChnlGrid *grn_uiz = (Ctrl->e.active)? (realChnlGrid *) calloc(Ctrl->nr, sizeof(*grn_uiz)) : NULL; - realChnlGrid *grn_uir = (Ctrl->e.active)? (realChnlGrid *) calloc(Ctrl->nr, sizeof(*grn_uir)) : NULL; + size_t ntot = ndepsrc * ndeprcv; + size_t idone = 0; + const char *modelname = grt_get_basename(modelpath); - K_INTEG_PROCESS KPROC = {0}; - prepare_static_grn( - mod1d, - Ctrl->nr, Ctrl->rs, - Ctrl->L.Length, - Ctrl->L.FIM.active ? Ctrl->L.FIM.Length : 0.0, - Ctrl->L.SAFIM.active ? Ctrl->L.SAFIM.tol : 0.0, - Ctrl->L.kcut, - Ctrl->K.k0, Ctrl->K.keps, Ctrl->K.use_kmax_ref, - Ctrl->C.convmet, - &KPROC); + for(size_t is = 0; is < ndepsrc; ++is){ + for(size_t ir = 0; ir < ndeprcv; ++ir){ + real_t zs = depsrcs[is]; + real_t zr = deprcvs[ir]; - //============================================================================== - // 计算静态格林函数 - grt_integ_static_grn( - mod1d, Ctrl->nr, Ctrl->rs, &KPROC, - Ctrl->e.active, grn, grn_uiz, grn_uir, - Ctrl->S.s_statsdir - ); - //============================================================================== - - real_t src_va = mod1d->Va[mod1d->isrc]; - real_t src_vb = mod1d->Vb[mod1d->isrc]; - real_t src_rho = mod1d->Rho[mod1d->isrc]; - real_t rcv_va = mod1d->Va[mod1d->ircv]; - real_t rcv_vb = mod1d->Vb[mod1d->ircv]; - real_t rcv_rho = mod1d->Rho[mod1d->ircv]; - - - // ================================================================================== - // 将结果保存为 nc 格式 - // ================================================================================== - int ncid, north_dimid, east_dimid; - const int ndims = 2; - int dimids[ndims]; - int north_varid, east_varid; - intChnlGrid u_varids; - intChnlGrid uiz_varids; - intChnlGrid uir_varids; - - // 创建 NC 文件 - NC_CHECK(nc_create(Ctrl->O.s_outgrid, NC_CLOBBER, &ncid)); - - // 写入全局属性 - NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "depsrc", NC_REAL, 1, &mod1d->depsrc)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "deprcv", NC_REAL, 1, &mod1d->deprcv)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "src_va", NC_REAL, 1, &src_va)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "src_vb", NC_REAL, 1, &src_vb)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "src_rho", NC_REAL, 1, &src_rho)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "rcv_va", NC_REAL, 1, &rcv_va)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "rcv_vb", NC_REAL, 1, &rcv_vb)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (ncid, NC_GLOBAL, "rcv_rho", NC_REAL, 1, &rcv_rho)); - // 是否计算了位移偏导也直接写到全局属性 - { - int tmp = Ctrl->e.active; - NC_CHECK(nc_put_att_int(ncid, NC_GLOBAL, "calc_upar", NC_INT, 1, &tmp)); - } + MODEL1D *mod1d = NULL; + if((mod1d = grt_read_mod1d_from_file(modelpath, zs, zr, false)) == NULL){ + exit(EXIT_FAILURE); + } + grt_set_mod1d_boundary(mod1d, topbound, botbound); - // 定义维度 - NC_CHECK(nc_def_dim(ncid, "north", Ctrl->X.nnorth, &north_dimid)); - NC_CHECK(nc_def_dim(ncid, "east", Ctrl->Y.neast, &east_dimid)); - dimids[0] = north_dimid; - dimids[1] = east_dimid; + // 拷贝模板，避免 integ 改写 cvgmet/dk 等影响后续深度 + K_INTEG_PROCESS local_K = *Kproc; + real_t hs = GRT_MAX(fabs(zs - zr), GRT_MIN_DEPTH_GAP_SRC_RCV); + local_K.k0 = k0 * PI / hs; - // 定义维度数组 - NC_CHECK(nc_def_var(ncid, "north", NC_REAL, 1, &north_dimid, &north_varid)); - NC_CHECK(nc_def_var(ncid, "east", NC_REAL, 1, &east_dimid, &east_varid)); + memset(grn, 0, nr * sizeof(*grn)); + if(calc_upar){ + memset(grn_uiz, 0, nr * sizeof(*grn_uiz)); + memset(grn_uir, 0, nr * sizeof(*grn_uir)); + } - // 定义不同震源不同分量的格林函数数组 - GRT_LOOP_ChnlGrid(im, c){ - int modr = GRT_SRC_M_ORDERS[im]; - char *s_title = NULL; + grt_integ_static_grn( + mod1d, nr, rs, &local_K, + calc_upar, grn, grn_uiz, grn_uir, statsstr); + + lib->src_va[is] = mod1d->Va[mod1d->isrc]; + lib->src_vb[is] = mod1d->Vb[mod1d->isrc]; + lib->src_rho[is] = mod1d->Rho[mod1d->isrc]; + lib->rcv_va[ir] = mod1d->Va[mod1d->ircv]; + lib->rcv_vb[ir] = mod1d->Vb[mod1d->ircv]; + lib->rcv_rho[ir] = mod1d->Rho[mod1d->ircv]; - if(modr==0 && GRT_ZRT_CODES[c]=='T') continue; + copy_grn_slice_with_sign(lib, is, ir, nr, calc_upar, grn, grn_uiz, grn_uir); - GRT_SAFE_ASPRINTF(&s_title, "%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); - NC_CHECK(nc_def_var(ncid, s_title, NC_REAL, ndims, dimids, &u_varids[im][c])); + idone++; + GRTRaiseInfo("[%zu/%zu] depsrc=%.6g deprcv=%.6g (%s) done.", + idone, ntot, zs, zr, modelname); - // 位移偏导 - if(Ctrl->e.active){ - GRT_SAFE_ASPRINTF(&s_title, "z%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); - NC_CHECK(nc_def_var(ncid, s_title, NC_REAL, ndims, dimids, &uiz_varids[im][c])); - GRT_SAFE_ASPRINTF(&s_title, "r%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); - NC_CHECK(nc_def_var(ncid, s_title, NC_REAL, ndims, dimids, &uir_varids[im][c])); + grt_free_mod1d(mod1d); } - GRT_SAFE_FREE_PTR(s_title); } - // 结束定义模式 - NC_CHECK(nc_enddef(ncid)); + // 将模型矩阵写入库，供后续 syn/应力按深度查层 + { + size_t nlayer = 0; + real_t (*modarr)[GRT_MODARR_NCOL] = grt_read_modarr_from_file(modelpath, &nlayer, false); + grt_stgrnlib_set_modarr(lib, nlayer, (const real_t (*)[GRT_MODARR_NCOL])modarr); + GRT_SAFE_FREE_PTR(modarr); + } - // 写入数据 - NC_CHECK(NC_FUNC_REAL(nc_put_var) (ncid, north_varid, Ctrl->X.norths)); - NC_CHECK(NC_FUNC_REAL(nc_put_var) (ncid, east_varid, Ctrl->Y.easts)); - real_t *tmpdata = (real_t *)calloc(Ctrl->nr, sizeof(real_t)); - GRT_LOOP_ChnlGrid(im, c){ - int modr = GRT_SRC_M_ORDERS[im]; + grt_stgrnlib_save_nc(lib, outpath); + GRTRaiseInfo("Static Green's function library saved in \"%s\".", outpath); - if(modr==0 && GRT_ZRT_CODES[c]=='T') continue; + GRT_SAFE_FREE_PTR(grn); + GRT_SAFE_FREE_PTR(grn_uiz); + GRT_SAFE_FREE_PTR(grn_uir); + grt_stgrnlib_free(lib); +} - int sgn0 = 1; - sgn0 = (GRT_ZRT_CODES[c]=='Z')? -1 : 1; - for(size_t ir=0; ir < Ctrl->nr; ++ir){ - tmpdata[ir] = sgn0 * grn[ir][im][c]; - } - NC_CHECK(NC_FUNC_REAL(nc_put_var) (ncid, u_varids[im][c], tmpdata)); +/** 子模块主函数 */ +int static_greenfn_main(int argc, char **argv){ + GRT_MODULE_CTRL *Ctrl = calloc(1, sizeof(*Ctrl)); - // 位移偏导 - if(Ctrl->e.active){ - for(size_t ir=0; ir < Ctrl->nr; ++ir){ - tmpdata[ir] = (-1) * sgn0 * grn_uiz[ir][im][c]; // 这里多乘的(-1)是因为对z的偏导，z需反向 - } - NC_CHECK(NC_FUNC_REAL(nc_put_var) (ncid, uiz_varids[im][c], tmpdata)); - for(size_t ir=0; ir < Ctrl->nr; ++ir){ - tmpdata[ir] = sgn0 * grn_uir[ir][im][c]; // 这里多乘的(-1)是因为对z的偏导，z需反向 - } - NC_CHECK(NC_FUNC_REAL(nc_put_var) (ncid, uir_varids[im][c], tmpdata)); + getopt_from_command(Ctrl, argc, argv); + + bool multi_depth = (Ctrl->D.ndepsrc > 1) || (Ctrl->D.ndeprcv > 1); + if(Ctrl->S.active){ + if(multi_depth){ + GRTRaiseWarning("-S is ignored for multi-depth STGRNLIB computation."); + } else { + // 单深度：stgrtstats/__ + GRT_SAFE_ASPRINTF(&Ctrl->S.s_statsdir, "stgrtstats"); + GRTCheckMakeDir(Ctrl->S.s_statsdir); + GRT_SAFE_ASPRINTF( + &Ctrl->S.s_statsdir, "%s/%s_%g_%g", + Ctrl->S.s_statsdir, + grt_get_basename(Ctrl->M.s_modelpath), + Ctrl->D.depsrcs[0], Ctrl->D.deprcvs[0]); + GRTCheckMakeDir(Ctrl->S.s_statsdir); } } - GRT_SAFE_FREE_PTR(tmpdata); - - // 关闭文件 - NC_CHECK(nc_close(ncid)); + K_INTEG_PROCESS KPROC = {0}; + prepare_static_grn( + Ctrl->nr, Ctrl->rs, + Ctrl->L.Length, + Ctrl->L.FIM.active ? Ctrl->L.FIM.Length : 0.0, + Ctrl->L.SAFIM.active ? Ctrl->L.SAFIM.tol : 0.0, + Ctrl->L.kcut, + Ctrl->K.keps, Ctrl->K.use_kmax_ref, + Ctrl->C.convmet, + &KPROC); - // 释放内存 - GRT_SAFE_FREE_PTR(grn); - GRT_SAFE_FREE_PTR(grn_uiz); - GRT_SAFE_FREE_PTR(grn_uir); + compute_stgrnlib_to_nc( + Ctrl->M.s_modelpath, + Ctrl->D.ndepsrc, Ctrl->D.depsrcs, + Ctrl->D.ndeprcv, Ctrl->D.deprcvs, + Ctrl->X.nnorth, Ctrl->X.norths, + Ctrl->Y.neast, Ctrl->Y.easts, + Ctrl->K.k0, + &KPROC, + Ctrl->B.topbound, Ctrl->B.botbound, + Ctrl->e.active, + Ctrl->O.s_outgrid, + Ctrl->S.s_statsdir); free_Ctrl(Ctrl); return EXIT_SUCCESS; diff --git a/pygrt/C_extension/src/static/grt_static_rotation.c b/pygrt/C_extension/src/static/grt_static_rotation.c index 50461383..5bd53d06 100644 --- a/pygrt/C_extension/src/static/grt_static_rotation.c +++ b/pygrt/C_extension/src/static/grt_static_rotation.c @@ -23,14 +23,21 @@ static void free_Ctrl(GRT_MODULE_CTRL *Ctrl){ static void print_help(){ printf("\n" "[grt static rotation] %s\n\n", GRT_VERSION);printf( -" Conbine spatial derivatives of static displacements\n" -" into rotation tensor, and write into the same nc file. \n" +" Combine spatial derivatives of static displacements\n" +" into rotation tensor, and write into the same nc file.\n" +" Input must be a static syn NetCDF computed with -e.\n" +" Both grid (-X/-Y) and points (-Q) layouts are supported.\n" " For example, \"ZR\" in variable names means\n" " 0.5*(u_{z,r} - u_{r,z}).\n" "\n\n" "Usage:\n" "----------------------------------------------------------------\n" " grt static rotation \n" +"\n" +"Examples:\n" +"----------------------------------------------------------------\n" +" grt static syn -Gstgrn.nc -Su1e16 -e -Ostsyn.nc\n" +" grt static rotation stsyn.nc\n" "\n\n\n" ); } @@ -51,28 +58,25 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ /** 由静态位移偏导合成旋转张量 */ static void compute_rotation( - size_t nnorth, size_t neast, const real_t *norths, const real_t *easts, + size_t npts, const real_t *norths, const real_t *easts, real_t *const u[GRT_CHANNEL_NUM], real_t *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], real_t *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE) { const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES; - for(size_t ix=0; ix 0.5*(u_{z,r} - u_{r,z}) GRT_SAFE_ASPRINTF(&s_title, "rotation_%c%c", toupper(chs[c]), toupper(chs[c2])); - NC_CHECK(nc_def_var(in_ncid, s_title, NC_REAL, ndims, in_dimids, &out_varids[c2][c])); + NC_CHECK(nc_def_var(in_ncid, s_title, NC_REAL, out_ndims, out_dimids, &out_varids[c2][c])); } GRT_SAFE_FREE_PTR(s_title); } @@ -169,25 +204,22 @@ int static_rotation_main(int argc, char **argv){ // 结束定义模式 NC_CHECK(nc_enddef(in_ncid)); - // 总震中距数 - size_t nr = nnorth * neast; - - // 先读入内存， + // 先读入内存 real_t *u[GRT_CHANNEL_NUM]; real_t *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; // 计算结果 real_t *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; for(int c=0; c\n" +"\n" +"Examples:\n" +"----------------------------------------------------------------\n" +" grt static syn -Gstgrn.nc -Su1e16 -e -Ostsyn.nc\n" +" grt static strain stsyn.nc\n" "\n\n\n" ); } @@ -49,32 +56,29 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ /** 由静态位移偏导合成应变张量 */ static void compute_strain( - size_t nnorth, size_t neast, const real_t *norths, const real_t *easts, + size_t npts, const real_t *norths, const real_t *easts, real_t *const u[GRT_CHANNEL_NUM], real_t *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], real_t *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], bool rot2ZNE) { const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES; - for(size_t ix=0; ix 0.5*(u_{z,r} + u_{r,z}) GRT_SAFE_ASPRINTF(&s_title, "strain_%c%c", toupper(chs[c]), toupper(chs[c2])); - NC_CHECK(nc_def_var(in_ncid, s_title, NC_REAL, ndims, in_dimids, &out_varids[c2][c])); + NC_CHECK(nc_def_var(in_ncid, s_title, NC_REAL, out_ndims, out_dimids, &out_varids[c2][c])); } GRT_SAFE_FREE_PTR(s_title); } @@ -171,25 +206,22 @@ int static_strain_main(int argc, char **argv){ // 结束定义模式 NC_CHECK(nc_enddef(in_ncid)); - // 总震中距数 - size_t nr = nnorth * neast; - - // 先读入内存， + // 先读入内存 real_t *u[GRT_CHANNEL_NUM]; real_t *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; // 计算结果 real_t *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; for(int c=0; c\n" +"\n" +"Examples:\n" +"----------------------------------------------------------------\n" +" grt static syn -Gstgrn.nc -Su1e16 -e -Ostsyn.nc\n" +" grt static stress stsyn.nc\n" "\n\n\n" ); } @@ -50,44 +59,59 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ /** 由静态位移偏导合成应力张量 */ static void compute_stress( - size_t nnorth, size_t neast, const real_t *norths, const real_t *easts, + size_t npts, const real_t *norths, const real_t *easts, real_t *const u[GRT_CHANNEL_NUM], real_t *const upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], real_t *const res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], - bool rot2ZNE, real_t mu, real_t lam) + bool rot2ZNE, const real_t *mu, const real_t *lam) { const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES; - for(size_t ix=0; ix mu*(u_{z,r} + u_{r,z}) GRT_SAFE_ASPRINTF(&s_title, "stress_%c%c", toupper(chs[c]), toupper(chs[c2])); - NC_CHECK(nc_def_var(in_ncid, s_title, NC_REAL, ndims, in_dimids, &out_varids[c2][c])); + NC_CHECK(nc_def_var(in_ncid, s_title, NC_REAL, out_ndims, out_dimids, &out_varids[c2][c])); } GRT_SAFE_FREE_PTR(s_title); } @@ -189,25 +263,22 @@ int static_stress_main(int argc, char **argv){ // 结束定义模式 NC_CHECK(nc_enddef(in_ncid)); - // 总震中距数 - size_t nr = nnorth * neast; - - // 先读入内存， + // 先读入内存 real_t *u[GRT_CHANNEL_NUM]; real_t *upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; // 计算结果 real_t *res[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; for(int c=0; cG.s_ingrid); + // C + grt_finite_fault_free(Ctrl->C.faults); + Ctrl->C.faults = NULL; + // X GRT_SAFE_FREE_PTR(Ctrl->X.norths); // Y GRT_SAFE_FREE_PTR(Ctrl->Y.easts); + // Q + GRT_SAFE_FREE_PTR(Ctrl->Q.s_path); + // O GRT_SAFE_FREE_PTR(Ctrl->O.s_outgrid); @@ -106,26 +136,57 @@ printf("\n" " + Transverse Clockwise (T),\n" " and the units are cm. You can add -N to rotate ZRT to ZNE.\n" "\n" -" You can also set -X/-Y to define a new north/east receiver grid;\n" -" synthesis interpolates Green's functions in epicentral distance.\n" -" (Green's functions themselves depend only on distance, even if the\n" -" input nc stores north/east dimensions for compatibility.)\n" +" Receivers: by default reuse the library north/east grid (from\n" +" greenfn -X/-Y or -R). Optionally redefine with -X/-Y (uniform\n" +" depth via -Dr when needed), or -Q for arbitrary points\n" +" (each with its own depth). -Q is mutually exclusive with -X/-Y\n" +" and -Dr. If the library was built with -R, the default grid is\n" +" a 1-D line (north=0, east=R); set -X/-Y or -Q to get a 2-D field.\n" +" Synthesis interpolates Green's functions in epicentral distance\n" +" (and in depth when needed).\n" "\n\n" "Usage:\n" "----------------------------------------------------------------\n" -" grt static syn -G -S[u] -O \n" +" # Point source\n" +" grt static syn -G -S[u] -O \n" +" [-Ds] [-Dr]\n" " [-M/[/]]\n" " [-T/////]\n" " [-F//] \n" -" [-X//] [-Y//]\n" +" [-X//] [-Y//] | [-Q]\n" " [-N] [-e] [-s]\n" "\n" +" # Finite faults (Coulomb format)\n" +" grt static syn -G -C[+i
/] -O\n" +" [-Dr] [-X//] [-Y//] | [-Q]\n" +" [-e] [-s]\n" +"\n" +" -G always points to a single 4D STGRNLIB nc file.\n" +" Depth options (without -Q) depend on the library shape:\n" +" ndepsrc=1, ndeprcv=1: do not set -Ds/-Dr; finite faults forbidden\n" +" ndepsrc=1, ndeprcv>1: -Dr required; do not set -Ds;\n" +" finite faults forbidden\n" +" ndepsrc>1, ndeprcv=1: -Ds required for point source;\n" +" do not set -Dr; finite faults allowed\n" +" ndepsrc>1, ndeprcv>1: -Ds required for point source;\n" +" -Dr required; finite faults allowed\n" +" With -Q, receiver depths come from the file; do not set -Dr.\n" +"\n" "\n\n" "Options:\n" "----------------------------------------------------------------\n" -" -G Filepath to input nc Green's Functions grid.\n" +" -G Filepath to a single STGRNLIB nc Green's function\n" +" library (dims depsrc×deprcv×north×east).\n" +"\n" +" -Ds Point-source source depth (km). Required when the\n" +" library has multiple source depths. Forbidden for\n" +" finite faults.\n" "\n" -" -S[u] Scale factor to all kinds of source. \n" +" -Dr Receiver depth (km) for grid receivers. Required\n" +" only when the library has multiple receiver depths\n" +" and -Q is not used; forbidden otherwise and with -Q.\n" +"\n" +" -S[u] Scale factor to all kinds of point source. \n" " + For Explosion, Shear and Moment Tensor,\n" " unit of is dyne-cm. \n" " + For Single Force, unit of is dyne.\n" @@ -136,8 +197,10 @@ printf("\n" "\n" " -O Filepath to output nc grid.\n" "\n" -" For source type, you can only set at most one of\n" +" For point source, you can only set at most one of\n" " '-M', '-T' and '-F'. If none, an Explosion is used.\n" +" For finite faults, use '-C' instead (mutually exclusive with\n" +" point-source options '-S'/'-M'/'-F'/'-T').\n" "\n" " -M/[/]\n" " Three angles to define a shear fault. \n" @@ -153,23 +216,44 @@ printf("\n" " North, East and Vertical(Downward) Forces.\n" " Notice they will be scaled by .\n" "\n" +" -C[+i
/]\n" +" Finite faults in Coulomb input format.\n" +" : fault file (skip two header lines).\n" +" Optional +i
/: along-strike / along-dip\n" +" subfault size (km). If omitted, both default to\n" +" min(dr, dz) of the Green's function library\n" +" (epicentral-distance and source-depth sampling).\n" +" Automatically enables -N (ZNE output).\n" +" Each fault: dip in (0, 90], bot > top (km).\n" +" Receiver locations default to the library grid;\n" +" optional -X/-Y or -Q to redefine.\n" +" Requires a library with ndepsrc>1.\n" +"\n" " -X//\n" " Set the equidistant points in the north direction.\n" " : start coordinate (km).\n" " : end coordinate (km).\n" " : sampling interval (km).\n" +" Mutually exclusive with -Q.\n" "\n" " -Y//\n" " Set the equidistant points in the east direction.\n" " : start coordinate (km).\n" " : end coordinate (km).\n" " : sampling interval (km).\n" +" Mutually exclusive with -Q.\n" +"\n" +" -Q Arbitrary receiver points from an ASCII file.\n" +" Each line: north east depth (km); lines starting\n" +" with # are comments. Mutually exclusive with\n" +" -X/-Y and -Dr (depths come from the file).\n" "\n" " -N Components of results will be Z, N, E.\n" "\n" -" -e Compute the spatial derivatives, ui_z and ui_r,\n" -" of displacement u. In filenames, prefix \"r\" means \n" -" ui_r and \"z\" means ui_z. \n" +" -e Also synthesize spatial derivatives of displacement.\n" +" Written as nc variables with prefixes z/r/t (ZRT)\n" +" or z/n/e (ZNE), e.g. zZ, rR. Required later for\n" +" static strain / stress / rotation.\n" "\n" " -s Silence all outputs.\n" "\n" @@ -177,26 +261,34 @@ printf("\n" "\n\n" "Examples:\n" "----------------------------------------------------------------\n" -" Say you have computed Static Green's functions with following command:\n" +" 2-D north/east grid (syn may omit -X/-Y and reuse the grid):\n" " grt static greenfn -Mmilrow -D2/0 -X-5/5/1 -Y-5/5/1 -Ostgrn.nc\n" -"\n" -" Then you can get static displacement of Explosion\n" " grt static syn -Gstgrn.nc -Su1e16 -Ostsyn_ex.nc\n" "\n" -" or Shear\n" +" Epicentral distances (recommended GF library). Syn should set\n" +" receivers with -X/-Y or -Q; otherwise output is a line along east:\n" +" grt static greenfn -Mmilrow -D2/0 -R0/7/0.1 -Ostgrn.nc\n" +" grt static syn -Gstgrn.nc -Su1e16 -X-5/5/0.5 -Y-5/5/0.5 -Ostsyn_ex.nc\n" +"\n" +" Other point sources (same -G file):\n" " grt static syn -Gstgrn.nc -Su1e16 -M100/20/80 -Ostsyn_dc.nc\n" +" grt static syn -Gstgrn.nc -Su1e16 -M100/20 -Ostsyn_ts.nc\n" +" grt static syn -Gstgrn.nc -S1e20 -F0.5/-1.2/3.3 -Ostsyn_sf.nc\n" +" grt static syn -Gstgrn.nc -Su1e16 -T2.3/0.2/-4.0/0.3/0.5/1.2 -Ostsyn_mt.nc\n" "\n" -" or Tension\n" -" grt static syn -Gstgrn.nc -Su1e16 -M100/20 -Ostsyn_dc.nc\n" +" Arbitrary receiver points (each line: north east depth in km):\n" +" grt static syn -Gstgrn.nc -Su1e16 -Qrcv.txt -Ostsyn_q.nc\n" "\n" -" or Single Force\n" -" grt static syn -Gstgrn.nc -S1e20 -F0.5/-1.2/3.3 -Ostsyn_sf.nc\n" +" Multi-depth library and interpolated source depth:\n" +" grt static greenfn -Mmilrow -Ds1,2,3 -Dr0 -R0/7/0.1 -Ostgrn.nc\n" +" grt static syn -Gstgrn.nc -Su1e16 -Ds1.5 -X-5/5/0.5 -Y-5/5/0.5 -Ostsyn.nc\n" "\n" -" or Moment Tensor\n" -" grt static syn -Gstgrn.nc -Su1e16 -T2.3/0.2/-4.0/0.3/0.5/1.2 -Ostsyn_mt.nc\n" +" Finite faults (Coulomb format; library must have ndepsrc>1):\n" +" grt static syn -Gstgrn.nc -Cfaults.inp+i1/1 -X-5/5/0.5 -Y-5/5/0.5 -Ostsyn_ff.nc\n" "\n" -" You can also set a new north/east receiver grid, for example\n" -" grt static syn -Gstgrn.nc -Su1e16 -Ostsyn_ex.nc -X-5/5/0.5 -Y-5/5/0.5\n" +" Spatial derivatives for later strain/stress/rotation:\n" +" grt static syn -Gstgrn.nc -Su1e16 -e -N -X-5/5/0.5 -Y-5/5/0.5 -Ostsyn.nc\n" +" grt static strain stsyn.nc\n" "\n\n\n" "\n" ); @@ -210,7 +302,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ sprintf(Ctrl->s_computeType, "%s", "EX"); int opt; - while ((opt = getopt(argc, argv, ":G:O:S:M:F:T:X:Y:Nesh")) != -1) { + while ((opt = getopt(argc, argv, ":G:O:S:M:F:T:C:X:Y:D:Q:Nesh")) != -1) { switch (opt) { // 输入 nc 文件名 case 'G': @@ -218,6 +310,29 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ Ctrl->G.s_ingrid = strdup(optarg); break; + // -Ds 或 -Dr（单值） + case 'D': + if(optarg[0] == 's'){ + Ctrl->D.s_active = true; + if(1 != sscanf(optarg + 1, "%lf", &Ctrl->D.depsrc)){ + GRTBadOptionError(Ds, ""); + } + if(Ctrl->D.depsrc < 0.0){ + GRTBadOptionError(Ds, "Negative source depth is not supported."); + } + } else if(optarg[0] == 'r'){ + Ctrl->D.r_active = true; + if(1 != sscanf(optarg + 1, "%lf", &Ctrl->D.deprcv)){ + GRTBadOptionError(Dr, ""); + } + if(Ctrl->D.deprcv < 0.0){ + GRTBadOptionError(Dr, "Negative receiver depth is not supported."); + } + } else { + GRTBadOptionError(D, "use -Ds or -Dr."); + } + break; + // 输出 nc 文件名 case 'O': Ctrl->O.active = true; @@ -308,6 +423,50 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } break; + // 从文件中读取有限断层 （Coulomb程序所用格式） -C[+i
/] + case 'C': + Ctrl->C.active = true; + { + Ctrl->computeType = GRT_SYN_DC; + sprintf(Ctrl->s_computeType, "%s", "DC"); + char *optarg_copy = strdup(optarg); + char *filepath = strtok(optarg_copy, "+"); + char *token = strtok(NULL, "+"); + if(filepath == NULL){ + GRT_SAFE_FREE_PTR(optarg_copy); + GRTBadOptionError(C, ""); + } + + // +i 可省略：dL/dW <= 0 表示按格林函数库步长取默认 + Ctrl->C.dL = 0.0; + Ctrl->C.dW = 0.0; + if(token != NULL){ + if(token[0] != 'i'){ + GRT_SAFE_FREE_PTR(optarg_copy); + GRTBadOptionError(C, ""); + } + if(2 != sscanf(token+1, "%lf/%lf", &Ctrl->C.dL, &Ctrl->C.dW)){ + GRT_SAFE_FREE_PTR(optarg_copy); + GRTBadOptionError(C, ""); + } + if(Ctrl->C.dL <= 0.0){ + GRT_SAFE_FREE_PTR(optarg_copy); + GRTBadOptionError(C, "dL(%f) <= 0.0", Ctrl->C.dL); + } + if(Ctrl->C.dW <= 0.0){ + GRT_SAFE_FREE_PTR(optarg_copy); + GRTBadOptionError(C, "dW(%f) <= 0.0", Ctrl->C.dW); + } + } + + char *fpath = strdup(filepath); + GRT_SAFE_FREE_PTR(optarg_copy); + + Ctrl->C.faults = grt_finite_fault_load_coulomb(fpath, &Ctrl->C.nfault); + GRT_SAFE_FREE_PTR(fpath); + } + break; + // X坐标数组，-Xx1/x2/dx case 'X': Ctrl->X.active = true; @@ -354,6 +513,12 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } break; + // 任意接收点文件，-Q + case 'Q': + Ctrl->Q.active = true; + Ctrl->Q.s_path = strdup(optarg); + break; + // 是否计算位移空间导数, 影响 calcUTypes 变量 case 'e': Ctrl->e.active = true; @@ -373,28 +538,51 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } } + // 要么通过 -S, -M/-F/-T 来指定点源，要么通过 -C 来指定有限断层 + bool isPointSource = Ctrl->S.active || Ctrl->M.active || Ctrl->F.active || Ctrl->T.active; + bool isFiniteFault = Ctrl->C.active; + if(isPointSource == isFiniteFault){ + GRTRaiseError("You may set either a point source or finite faults — not both, and not neither. Use \"-h\" for help.\n"); + } + // 检查必选项有没有设置 GRTCheckOptionSet(argc > 1); GRTCheckOptionActive(Ctrl, G); GRTCheckOptionActive(Ctrl, O); - GRTCheckOptionActive(Ctrl, S); + if(isPointSource) GRTCheckOptionActive(Ctrl, S); - // 只能使用一种震源 - if(Ctrl->M.active + Ctrl->F.active + Ctrl->T.active > 1){ + // 点源只能使用一种震源 + if(isPointSource && (Ctrl->M.active + Ctrl->F.active + Ctrl->T.active > 1)){ GRTRaiseError("Only support at most one of \"-M\", \"-F\" and \"-T\". Use \"-h\" for help.\n"); } - // 指定新网格时必须同时指定 -X 和 -Y + // -Q 与 -X/-Y 互斥 + if(Ctrl->Q.active && (Ctrl->X.active || Ctrl->Y.active)){ + GRTRaiseError("\"-Q\" is mutually exclusive with \"-X\"/\"-Y\". Use \"-h\" for help.\n"); + } + + // 指定新接收点网格时必须同时指定 -X 和 -Y if(Ctrl->X.active ^ Ctrl->Y.active){ GRTRaiseError("If you want to set a new north/east grid, you need set \"-X\" and \"-Y\" both.\n"); } Ctrl->isnewNEgrid = Ctrl->X.active; + + // -Q 时深度来自文件，禁止 -Dr + if(Ctrl->Q.active && Ctrl->D.r_active){ + GRTRaiseError("Do not set -Dr with -Q; receiver depths come from the points file.\n"); + } + + // 有限断层：自动启用 -N；接收点默认延用库坐标，可用 -X/-Y 或 -Q 覆盖 + if(isFiniteFault){ + Ctrl->N.active = true; + } + + Ctrl->isPointSource = isPointSource; + Ctrl->isFiniteFault = isFiniteFault; } -/** - * 在单个震中距点上，由静态格林函数合成三分量（及可选空间偏导） - */ +/** 在单个震中距点上，由静态格林函数合成三分量（及可选空间偏导） */ static void static_syn_from_gf_one( real_t azrad, size_t ir_pick, real_t dist0, const realChnlGrid *u, const realChnlGrid *uiz, const realChnlGrid *uir, @@ -469,354 +657,743 @@ static void static_syn_from_gf_one( /** * 由静态格林函数合成三分量位移场（及可选空间偏导） * - * 输入 nc 虽以 north/east 存储，合成时按 r=hypot(north,east) 做一维震中距插值； - * 可换到新的接收点 north/east 网格。r=0 时强制方位角为 0（e_r→N、e_θ→E） + * GF 侧使用已准备好的升序震中距元数据 + * （sort_rs0 / sort_rs0_idx / isUniform / dr）；查询点为平坦 north/east 列表 + * r=0 时强制方位角为 0（e_r→N、e_θ→E） * * 数组布局：u[采样点][震源][分量]、syn[接收点][分量]、 * syn_upar[接收点][偏导方向][分量]。uiz/uir 在 calc_upar=false 时可传 NULL */ static void static_syn_from_gf( - size_t nnorth0, const real_t *norths0, size_t neast0, const real_t *easts0, - size_t nnorth, const real_t *norths, size_t neast, const real_t *easts, + size_t nr0, const real_t *sort_rs0, const size_t *sort_rs0_idx, + bool isUniform, real_t dr, + size_t npts, const real_t *norths, const real_t *easts, const realChnlGrid *u, const realChnlGrid *uiz, const realChnlGrid *uir, GRT_SYN_TYPE computeType, real_t M0, real_t VpVs_ratio, const real_t mchn[GRT_MECHANISM_NUM], bool rot2ZNE, bool calc_upar, real_t (*syn)[GRT_CHANNEL_NUM], real_t (*syn_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]) { - size_t nr0 = nnorth0 * neast0; - - // 由 north/east 派生震中距序列（GF 只依赖 r） - real_t *rs0 = (real_t *)calloc(nr0, sizeof(real_t)); - real_t *sort_rs0 = (real_t *)calloc(nr0, sizeof(real_t)); - size_t *sort_rs0_idx = (size_t *)calloc(nr0, sizeof(size_t)); - for(size_t inorth = 0; inorth < nnorth0; ++inorth){ - for(size_t ieast = 0; ieast < neast0; ++ieast){ - size_t idx = ieast + inorth*neast0; - rs0[idx] = hypot(norths0[inorth], easts0[ieast]); - sort_rs0_idx[idx] = idx; - } - } - memcpy(sort_rs0, rs0, nr0 * sizeof(*sort_rs0)); + // 每个接收点逐个处理 + for(size_t ir = 0; ir < npts; ++ir){ + real_t north = norths[ir]; + real_t east = easts[ir]; + + real_t dist = hypot(north, east); + + // 方位角；r=0 时 atan2(0,0) 无定义，约定 e_r→N、e_θ→E ⇒ az=0 + real_t azrad = GRT_IS_ZERO(dist) ? 0.0 : atan2(east, north); + + // syn/syn_upar 不在此处清零，以便多次调用（有限断层子源）时对 syn 做累加 + // 调用方需保证首次调用前 syn/syn_upar 已清零（如 calloc） - // 还未排序前，先判断是否是一个等距升序数组，这样对于加快后续查找 - bool isUniform = (nr0 > 2); - real_t dr = (nr0 > 1)? rs0[1] - rs0[0] : 0.0; - for(size_t ir = 1; ir < nr0-1; ++ir){ - if( fabs(2.0*rs0[ir] - (rs0[ir-1] + rs0[ir+1])) > 1e-3 || rs0[ir-1] >= rs0[ir] || rs0[ir] >= rs0[ir+1]){ - isUniform = false; - break; + // 检查是否越界（允许查询点为精确的零震中距） + bool r_OutofBound = (dist < sort_rs0[0] - 1e-8 || dist > sort_rs0[nr0-1] + 1e-8); + if(r_OutofBound){ + GRTRaiseWarning("(north, east)=(%.3e, %.3e) is out of distance bounds, skip.", north, east); + continue; } - } - if(! isUniform){ - if(grt_argsort(rs0, nr0, sizeof(*rs0), grt_compare_real_t, sort_rs0_idx) != 0){ - GRTRaiseError("Unable to sort source-grid distances."); + size_t sort_ir_pick = 0, sort_ir_pick1 = 0; + if(isUniform){ + sort_ir_pick = (size_t)((dist - sort_rs0[0]) / dr); + } else { + for(sort_ir_pick = 0; sort_ir_pick < nr0-1; ++sort_ir_pick) if(sort_rs0[sort_ir_pick+1] > dist) break; } + sort_ir_pick1 = GRT_MIN(sort_ir_pick + 1, nr0-1); + + // 重复震中距时避免除零（-X/-Y 建库时对称点可能 r 相同） + real_t r0 = sort_rs0[sort_ir_pick]; + real_t r1 = sort_rs0[sort_ir_pick1]; + real_t drs = (sort_ir_pick == sort_ir_pick1 || fabs(r1 - r0) < 1e-15) + ? 0.0 : (dist - r0) / (r1 - r0); - for(size_t i = 0; i < nr0; ++i){ - sort_rs0[i] = rs0[sort_rs0_idx[i]]; + real_t syn2[GRT_CHANNEL_NUM] = {0.0}, syn2_upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM] = {{0.0}}; + + size_t iir[2] = {sort_ir_pick, sort_ir_pick1}; + real_t facr[2] = {1.0 - drs, drs}; + for(int j = 0; j < 2; ++j){ + if(j==1 && (sort_ir_pick == sort_ir_pick1 || fabs(r1 - r0) < 1e-15)) continue; + + size_t ir_pick = sort_rs0_idx[iir[j]]; + real_t dist0 = sort_rs0[iir[j]]; + static_syn_from_gf_one(azrad, ir_pick, dist0, u, uiz, uir, computeType, M0, VpVs_ratio, mchn, rot2ZNE, calc_upar, syn2, syn2_upar); + + for(int c = 0; c < GRT_CHANNEL_NUM; ++c){ + syn[ir][c] += facr[j] * syn2[c]; + for(int c2 = 0; c2 < GRT_CHANNEL_NUM; ++c2){ + syn_upar[ir][c][c2] += facr[j] * syn2_upar[c][c2]; + } + } } } +} - // 每个接收点逐个处理 - for(size_t inorth = 0; inorth < nnorth; ++inorth){ - real_t north = norths[inorth]; - for(size_t ieast = 0; ieast < neast; ++ieast){ - real_t east = easts[ieast]; - // 震中距 - real_t dist = hypot(north, east); +/** + * depsrc×deprcv 双线性角点循环：对给定水平坐标子集做震中距合成，按 fac 累加到 syn + * + * 角点权 fac = w_src * w_rcv；结果写入 syn[ipt0 + i]（i = 0..nloc-1） + * + * @param[in] lib 静态格林函数库 + * @param[in] is0 depsrc 括号左端下标 + * @param[in] is1 depsrc 括号右端下标 + * @param[in] ws depsrc 插值权（落在 is1 侧） + * @param[in] na depsrc 角点个数（1 或 2） + * @param[in] ir0 deprcv 括号左端下标 + * @param[in] ir1 deprcv 括号右端下标 + * @param[in] wr deprcv 插值权（落在 ir1 侧） + * @param[in] ipt0 写回 syn / syn_upar 的起始下标 + * @param[in] nloc 本次合成的接收点数 + * @param[in] loc_n 水平 north 坐标 (km)，长度 nloc + * @param[in] loc_e 水平 east 坐标 (km)，长度 nloc + * @param[in] computeType 震源类型 + * @param[in] M0 标量矩或 potency（见 scale_by_src_mu） + * @param[in] scale_by_src_mu 为真时用角点 μ 将 M0 转为矩 + * @param[in] mchn 震源机制参数数组 + * @param[in] rot2ZNE 是否输出 ZNE + * @param[in] calc_upar 是否合成位移空间偏导 + * @param[out] tmp 单角点位移缓冲，容量 >= nloc + * @param[out] tmp_upar 单角点偏导缓冲，容量 >= nloc + * @param[in,out] syn 位移累加输出，下标从 ipt0 起 + * @param[in,out] syn_upar 偏导累加输出，下标从 ipt0 起；仅在 calc_upar 为真时写入 + */ +static void static_syn_ps_depth_corners( + const STGRNLIB *lib, + size_t is0, size_t is1, real_t ws, int na, + size_t ir0, size_t ir1, real_t wr, + size_t ipt0, size_t nloc, const real_t *loc_n, const real_t *loc_e, + GRT_SYN_TYPE computeType, real_t M0, bool scale_by_src_mu, + const real_t mchn[GRT_MECHANISM_NUM], + bool rot2ZNE, bool calc_upar, + real_t (*tmp)[GRT_CHANNEL_NUM], + real_t (*tmp_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], + real_t (*syn)[GRT_CHANNEL_NUM], + real_t (*syn_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]) +{ + size_t is_idx[2] = {is0, is1}; + size_t ir_idx[2] = {ir0, ir1}; + int nb = (ir0 == ir1) ? 1 : 2; + + for(int a = 0; a < na; ++a){ + for(int b = 0; b < nb; ++b){ + real_t fac = ((a == 0) ? (1.0 - ws) : ws) * ((b == 0) ? (1.0 - wr) : wr); + if(fabs(fac) < 1e-15) continue; + + memset(tmp, 0, nloc * sizeof(real_t) * GRT_CHANNEL_NUM); + if(calc_upar){ + memset(tmp_upar, 0, nloc * sizeof(real_t) * GRT_CHANNEL_NUM * GRT_CHANNEL_NUM); + } - // 方位角；r=0 时 atan2(0,0) 无定义，约定 e_r→N、e_θ→E ⇒ az=0 - real_t azrad = GRT_IS_ZERO(dist) ? 0.0 : atan2(east, north); + size_t is = is_idx[a]; + size_t ir = ir_idx[b]; + + // 角点介质取 depsrcs[is] + real_t va = lib->src_va[is]; + real_t vb = lib->src_vb[is]; + real_t rho = lib->src_rho[is]; + real_t VpVs_ratio = (vb == 0.0) ? 0.0 : (va / vb); + real_t M0_use = M0; + if(scale_by_src_mu){ + M0_use = M0 * (vb * vb * rho * 1e10); // dyne/cm^2 * potency + } - size_t ir = ieast + inorth * neast; + static_syn_from_gf( + lib->nr, lib->sort_rs, lib->sort_rs_idx, lib->isUniform, lib->dr, + nloc, loc_n, loc_e, + lib->u[is][ir], + calc_upar ? lib->uiz[is][ir] : NULL, + calc_upar ? lib->uir[is][ir] : NULL, + computeType, M0_use, VpVs_ratio, mchn, + rot2ZNE, calc_upar, + tmp, tmp_upar + ); + + for(size_t i = 0; i < nloc; ++i){ + size_t ipt = ipt0 + i; + for(int c = 0; c < GRT_CHANNEL_NUM; ++c){ + syn[ipt][c] += fac * tmp[i][c]; + if(calc_upar){ + for(int c2 = 0; c2 < GRT_CHANNEL_NUM; ++c2){ + syn_upar[ipt][c][c2] += fac * tmp_upar[i][c][c2]; + } + } + } + } + } + } +} - memset(syn[ir], 0, sizeof(syn[ir])); - memset(syn_upar[ir], 0, sizeof(syn_upar[ir])); - // 检查是否越界（允许查询点为精确的零震中距） - bool r_OutofBound = (dist < sort_rs0[0] - 1e-8 || dist > sort_rs0[nr0-1] + 1e-8); - if(r_OutofBound){ - GRTRaiseWarning("(north, east)=(%.3e, %.3e) is out of distance bounds, skip.", north, east); - continue; - } +/** + * 基于 STGRNLIB：对各 depsrc×deprcv 邻点做震中距 1D 合成，再对结果做 2D 组合 + * + * shared_depth 为真（-X/-Y 网格或延用库水平网格）：全部接收点共面，用 depths[0] 求一次 deprcv 括号后批量合成 + * shared_depth 为假（-Q 任意点）：逐点求 deprcv 括号并合成，不做深度归组 + * 各角点的 Vp/Vs（及可选 μ）取自对应 depsrcs 采样； + * scale_by_src_mu 为真时，M0 为 potency，角点矩为 M0 * μ[is] + * + * 输出 syn / syn_upar 按原 npts 下标累加，调用方需事先清零（如 calloc） + */ +static void static_syn_from_gf_PS( + const STGRNLIB *lib, real_t depsrc, + size_t npts, const real_t *norths, const real_t *easts, const real_t *depths, + bool shared_depth, + GRT_SYN_TYPE computeType, real_t M0, bool scale_by_src_mu, + const real_t mchn[GRT_MECHANISM_NUM], + bool rot2ZNE, bool calc_upar, + real_t (*syn)[GRT_CHANNEL_NUM], + real_t (*syn_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]) +{ + if(lib == NULL || lib->ndepsrc == 0 || lib->ndeprcv == 0){ + GRTRaiseError("empty STGRNLIB."); + } + if(npts == 0 || norths == NULL || easts == NULL || depths == NULL){ + GRTRaiseError("empty receiver points."); + } + if(calc_upar && !lib->calc_upar){ + GRTRaiseError("STGRNLIB has no displacement derivatives, cannot set calc_upar."); + } - size_t sort_ir_pick = 0, sort_ir_pick1 = 0; - if(isUniform){ - sort_ir_pick = (size_t)((dist - sort_rs0[0]) / dr); - } else { - for(sort_ir_pick = 0; sort_ir_pick < nr0-1; ++sort_ir_pick) if(sort_rs0[sort_ir_pick+1] > dist) break; + // 震源深度括号：depsrc 落在 depsrcs[is0], depsrcs[is1] 之间，权为 ws + size_t is0, is1; + real_t ws; + if(!grt_locateLinearInterp(lib->depsrcs, lib->ndepsrc, depsrc, &is0, &is1, &ws)){ + GRTRaiseError( + "Source depth %.6g km is out of Green's function depsrc range [%.6g, %.6g].", + depsrc, lib->depsrcs[0], lib->depsrcs[lib->ndepsrc - 1]); + } + int na = (is0 == is1) ? 1 : 2; // 落在采样点上时只取一侧 + + // 单角点合成缓冲：共面时容量 npts，逐点时容量 1 + size_t nbuf = shared_depth ? npts : 1; + real_t (*tmp)[GRT_CHANNEL_NUM] = + (real_t (*)[GRT_CHANNEL_NUM])calloc(nbuf, sizeof(real_t) * GRT_CHANNEL_NUM); + real_t (*tmp_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM] = + (real_t (*)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM])calloc( + nbuf, sizeof(real_t) * GRT_CHANNEL_NUM * GRT_CHANNEL_NUM); + + if(shared_depth){ + // 网格：统一深度 depths[0]，一次括号 + 批量震中距合成 + real_t zr = depths[0]; + size_t ir0, ir1; + real_t wr; + if(!grt_locateLinearInterp(lib->deprcvs, lib->ndeprcv, zr, &ir0, &ir1, &wr)){ + GRTRaiseError( + "Receiver depth %.6g km is out of Green's function deprcv range [%.6g, %.6g].", + zr, lib->deprcvs[0], lib->deprcvs[lib->ndeprcv - 1]); + } + static_syn_ps_depth_corners( + lib, is0, is1, ws, na, ir0, ir1, wr, + 0, npts, norths, easts, + computeType, M0, scale_by_src_mu, mchn, + rot2ZNE, calc_upar, tmp, tmp_upar, syn, syn_upar); + } else { + // 任意点：逐点括号与合成 + for(size_t ipt = 0; ipt < npts; ++ipt){ + real_t zr = depths[ipt]; + size_t ir0, ir1; + real_t wr; + if(!grt_locateLinearInterp(lib->deprcvs, lib->ndeprcv, zr, &ir0, &ir1, &wr)){ + GRTRaiseError( + "Receiver depth %.6g km is out of Green's function deprcv range [%.6g, %.6g].", + zr, lib->deprcvs[0], lib->deprcvs[lib->ndeprcv - 1]); } - sort_ir_pick1 = GRT_MIN(sort_ir_pick + 1, nr0-1); + static_syn_ps_depth_corners( + lib, is0, is1, ws, na, ir0, ir1, wr, + ipt, 1, &norths[ipt], &easts[ipt], + computeType, M0, scale_by_src_mu, mchn, + rot2ZNE, calc_upar, tmp, tmp_upar, syn, syn_upar); + } + } - real_t r0 = sort_rs0[sort_ir_pick]; - real_t r1 = sort_rs0[sort_ir_pick1]; - // 重复震中距时避免除零（-X/-Y 建库时对称点可能 r 相同） - real_t drs = (sort_ir_pick == sort_ir_pick1 || fabs(r1 - r0) < 1e-15) - ? 0.0 : (dist - r0) / (r1 - r0); + GRT_SAFE_FREE_PTR(tmp); + GRT_SAFE_FREE_PTR(tmp_upar); +} - real_t syn2[GRT_CHANNEL_NUM] = {0.0}, syn2_upar[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM] = {{0.0}}; - size_t iir[2] = {sort_ir_pick, sort_ir_pick1}; - real_t facr[2] = {1.0 - drs, drs}; - for(int j = 0; j < 2; ++j){ - if(j==1 && (sort_ir_pick == sort_ir_pick1 || fabs(r1 - r0) < 1e-15)) continue; +/** + * 单个有限断层：子源以 potency 传入，各 depsrc 邻点用对应 μ 合成后再 2D 组合 + */ +static void static_syn_one_finite_fault( + const STGRNLIB *lib, const FINITE_FAULT *fault, + real_t dL, real_t dW, real_t W, real_t L, size_t nW, size_t nL, + size_t npts, const real_t *norths, const real_t *easts, const real_t *depths, + bool shared_depth, bool calc_upar, + real_t (*syn)[GRT_CHANNEL_NUM], + real_t (*syn_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]) +{ + real_t mchn[GRT_MECHANISM_NUM] = {0}; + mchn[0] = fault->strike; + mchn[1] = fault->dip; + mchn[2] = fault->rake; + + // 按子源并行：各线程累加到私有缓冲，最后归约到 syn，避免对同一接收点写竞争 + size_t nsub = nW * nL; + #pragma omp parallel default(shared) if(nsub > 1) + { + real_t *rcv_norths = (real_t *)calloc(npts, sizeof(real_t)); + real_t *rcv_easts = (real_t *)calloc(npts, sizeof(real_t)); + real_t (*local_syn)[GRT_CHANNEL_NUM] = + (real_t (*)[GRT_CHANNEL_NUM])calloc(npts, sizeof(real_t) * GRT_CHANNEL_NUM); + real_t (*local_syn_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM] = + (real_t (*)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM])calloc( + npts, sizeof(real_t) * GRT_CHANNEL_NUM * GRT_CHANNEL_NUM); + + #pragma omp for collapse(2) schedule(guided) + for(size_t iW = 0; iW < nW; ++iW){ + for(size_t iL = 0; iL < nL; ++iL){ + FINITE_SUBFAULT sub; + grt_finite_fault_subfault(fault, dL, dW, W, L, iW, iL, &sub); + + for(size_t ipt = 0; ipt < npts; ++ipt){ + rcv_norths[ipt] = norths[ipt] - sub.north; + rcv_easts[ipt] = easts[ipt] - sub.east; + } - size_t ir_pick = sort_rs0_idx[iir[j]]; - real_t dist0 = sort_rs0[iir[j]]; - static_syn_from_gf_one(azrad, ir_pick, dist0, u, uiz, uir, computeType, M0, VpVs_ratio, mchn, rot2ZNE, calc_upar, syn2, syn2_upar); + static_syn_from_gf_PS( + lib, sub.depsrc, + npts, rcv_norths, rcv_easts, depths, + shared_depth, + GRT_SYN_DC, sub.potency, true, mchn, + true, calc_upar, + local_syn, local_syn_upar + ); + } + } + #pragma omp critical(static_syn_ff_reduce) + { + for(size_t ipt = 0; ipt < npts; ++ipt){ for(int c = 0; c < GRT_CHANNEL_NUM; ++c){ - syn[ir][c] += facr[j] * syn2[c]; - for(int c2 = 0; c2 < GRT_CHANNEL_NUM; ++c2){ - syn_upar[ir][c][c2] += facr[j] * syn2_upar[c][c2]; + syn[ipt][c] += local_syn[ipt][c]; + if(calc_upar){ + for(int c2 = 0; c2 < GRT_CHANNEL_NUM; ++c2){ + syn_upar[ipt][c][c2] += local_syn_upar[ipt][c][c2]; + } } } } } - } - GRT_SAFE_FREE_PTR(rs0); - GRT_SAFE_FREE_PTR(sort_rs0); - GRT_SAFE_FREE_PTR(sort_rs0_idx); + GRT_SAFE_FREE_PTR(rcv_norths); + GRT_SAFE_FREE_PTR(rcv_easts); + GRT_SAFE_FREE_PTR(local_syn); + GRT_SAFE_FREE_PTR(local_syn_upar); + } } -/** 子模块主函数 */ -int static_syn_main(int argc, char **argv){ - GRT_MODULE_CTRL *Ctrl = calloc(1, sizeof(*Ctrl)); +/** + * 有限断层合成（Coulomb 格式断层 + STGRNLIB） + * dL、dW 均 <=0 时取 grt_stgrnlib_default_subfault_size(lib) + * shared_depth 含义同 static_syn_from_gf_PS（网格共面 / 任意点逐点） + */ +static void static_syn_from_gf_FF( + const STGRNLIB *lib, + size_t nfault, const FINITE_FAULT *faults, + real_t dL, real_t dW, + size_t npts, const real_t *norths, const real_t *easts, const real_t *depths, + bool shared_depth, bool calc_upar, + real_t (*syn)[GRT_CHANNEL_NUM], + real_t (*syn_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]) +{ + if(lib == NULL || lib->ndepsrc == 0 || lib->ndeprcv == 0){ + GRTRaiseError("empty STGRNLIB."); + } + if(lib->ndepsrc <= 1){ + GRTRaiseError("Finite faults require a Green's function library with ndepsrc > 1."); + } + if(npts == 0 || norths == NULL || easts == NULL || depths == NULL){ + GRTRaiseError("empty receiver points."); + } + if((dL <= 0.0) != (dW <= 0.0)){ + GRTRaiseError("set both dL and dW, or neither for default."); + } + if(dL <= 0.0){ + dL = dW = grt_stgrnlib_default_subfault_size(lib); + GRTRaiseInfo("finite fault: use default dL = dW = %.6g km", dL); + } + if(calc_upar && !lib->calc_upar){ + GRTRaiseError("STGRNLIB has no displacement derivatives, cannot set calc_upar."); + } - getopt_from_command(Ctrl, argc, argv); + // 预先检查所有接收深度均在库范围内 + { + size_t i0, i1; + real_t w; + for(size_t ipt = 0; ipt < npts; ++ipt){ + if(!grt_locateLinearInterp(lib->deprcvs, lib->ndeprcv, depths[ipt], &i0, &i1, &w)){ + GRTRaiseError( + "Receiver depth %.6g km is out of Green's function deprcv range [%.6g, %.6g].", + depths[ipt], lib->deprcvs[0], lib->deprcvs[lib->ndeprcv - 1]); + } + } + } - // 输出分量格式，即是否需要旋转到ZNE - bool rot2ZNE = Ctrl->N.active; + for(size_t ifault = 0; ifault < nfault; ++ifault){ + FINITE_FAULT f = faults[ifault]; + grt_finite_fault_set_derived(&f); - // 根据参数设置，选择分量名 - const char *chs = (rot2ZNE)? GRT_ZNE_CODES : GRT_ZRT_CODES; + real_t W, L; + size_t nW, nL; + grt_finite_fault_subdiv(&f, dL, dW, &W, &L, &nW, &nL); + GRTRaiseInfo("finite fault[%zu/%zu]: nsubfaults = %zu", ifault + 1, nfault, nW * nL); - // nc 文件相关变量 - int in_ncid; - int in_x_dimid, in_y_dimid; - int in_x_varid, in_y_varid; - const int ndims = 2; - intChnlGrid in_u_varids; - intChnlGrid in_uiz_varids; - intChnlGrid in_uir_varids; - int out_ncid; - int out_x_dimid, out_y_dimid; - int out_x_varid, out_y_varid; - int out_dimids[ndims]; - int out_syn_varids[GRT_CHANNEL_NUM]; - int out_syn_upar_varids[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; - - size_t nnorth, neast, nnorth0, neast0; - real_t *norths, *easts, *norths0, *easts0; - size_t nr, nr0; - - // 打开 nc 文件 - GRTCheckFileExist(Ctrl->G.s_ingrid); - NC_CHECK(nc_open(Ctrl->G.s_ingrid, NC_NOWRITE, &in_ncid)); - NC_CHECK(nc_create(Ctrl->O.s_outgrid, NC_CLOBBER, &out_ncid)); - - // 读取全局属性，视情况计算 src_mu - real_t src_va=0.0, src_vb=0.0, src_rho=0.0, src_mu=0.0; - NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "src_va", &src_va)); - NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "src_vb", &src_vb)); - NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "src_rho", &src_rho)); - src_mu = src_vb*src_vb*src_rho*1e10; - if(Ctrl->S.mult_src_mu) Ctrl->S.M0 *= src_mu; - - // 读入的数据是否有位移偏导 - int calc_upar; - NC_CHECK(nc_get_att_int(in_ncid, NC_GLOBAL, "calc_upar", &calc_upar)); - if(Ctrl->e.active && calc_upar == 0){ - GRTRaiseError("Input grid didn't have displacement derivatives, you can't set -e."); + static_syn_one_finite_fault( + lib, &f, dL, dW, W, L, nW, nL, + npts, norths, easts, depths, + shared_depth, calc_upar, syn, syn_upar + ); } +} - // 复制属性 - NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "src_va", NC_REAL, 1, &src_va)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "src_vb", NC_REAL, 1, &src_vb)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "src_rho", NC_REAL, 1, &src_rho)); - NC_CHECK(nc_put_att_int(out_ncid, NC_GLOBAL, "calc_upar", NC_INT, 1, &calc_upar)); - { - real_t rcv_va=0.0, rcv_vb=0.0, rcv_rho=0.0; - NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "rcv_va", &rcv_va)); - NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "rcv_vb", &rcv_vb)); - NC_CHECK(NC_FUNC_REAL(nc_get_att) (in_ncid, NC_GLOBAL, "rcv_rho", &rcv_rho)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "rcv_va", NC_REAL, 1, &rcv_va)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "rcv_vb", NC_REAL, 1, &rcv_vb)); - NC_CHECK(NC_FUNC_REAL(nc_put_att) (out_ncid, NC_GLOBAL, "rcv_rho", NC_REAL, 1, &rcv_rho)); - } - // 是否旋转到ZNE记录到全局属性 - { - int rot2ZNE_int = rot2ZNE; - NC_CHECK(nc_put_att_int(out_ncid, NC_GLOBAL, "rot2ZNE", NC_INT, 1, &rot2ZNE_int)); - } - - // 震源类型写入全局属性 - NC_CHECK(nc_put_att_text(out_ncid, NC_GLOBAL, "computeType", strlen(Ctrl->s_computeType), Ctrl->s_computeType)); - - // 读入坐标变量 dimid, varid - NC_CHECK(nc_inq_dimid(in_ncid, "north", &in_x_dimid)); - NC_CHECK(nc_inq_dimlen(in_ncid, in_x_dimid, &nnorth0)); - NC_CHECK(nc_inq_dimid(in_ncid, "east", &in_y_dimid)); - NC_CHECK(nc_inq_dimlen(in_ncid, in_y_dimid, &neast0)); - norths0 = (real_t *)calloc(nnorth0, sizeof(real_t)); - easts0 = (real_t *)calloc(neast0, sizeof(real_t)); - - // 根据情况使用 -X/-Y 指定的新接收点网格 - nnorth = (Ctrl->isnewNEgrid)? Ctrl->X.nnorth : nnorth0; - neast = (Ctrl->isnewNEgrid)? Ctrl->Y.neast : neast0; - norths = (Ctrl->isnewNEgrid)? Ctrl->X.norths : norths0; - easts = (Ctrl->isnewNEgrid)? Ctrl->Y.easts : easts0; - - // 写入坐标变量 dimid, varid - NC_CHECK(nc_def_dim(out_ncid, "north", nnorth, &out_x_dimid)); - NC_CHECK(nc_def_dim(out_ncid, "east", neast, &out_y_dimid)); - NC_CHECK(nc_def_var(out_ncid, "north", NC_REAL, 1, &out_x_dimid, &out_x_varid)); - NC_CHECK(nc_def_var(out_ncid, "east", NC_REAL, 1, &out_y_dimid, &out_y_varid)); - out_dimids[0] = out_x_dimid; - out_dimids[1] = out_y_dimid; - - // 读入格林函数 varid - GRT_LOOP_ChnlGrid(im, c){ - int modr = GRT_SRC_M_ORDERS[im]; - char *s_title = NULL; - if(modr==0 && GRT_ZRT_CODES[c]=='T') continue; +/** 按库形态校验 -Ds/-Dr/-C（先禁止项，再必填项） */ +static void check_syn_depth_options(const GRT_MODULE_CTRL *Ctrl, const STGRNLIB *lib) +{ + bool multi_src = (lib->ndepsrc > 1); + bool multi_rcv = (lib->ndeprcv > 1); - GRT_SAFE_ASPRINTF(&s_title, "%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); - NC_CHECK(nc_inq_varid(in_ncid, s_title, &in_u_varids[im][c])); + if(Ctrl->isFiniteFault && !multi_src){ + GRTRaiseError("Finite faults require a Green's function library with ndepsrc > 1."); + } + if(Ctrl->isFiniteFault && Ctrl->D.s_active){ + GRTRaiseError("Do not set -Ds for finite faults; source depths come from the fault geometry."); + } - // 位移偏导 - if(Ctrl->e.active){ - GRT_SAFE_ASPRINTF(&s_title, "z%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); - NC_CHECK(nc_inq_varid(in_ncid, s_title, &in_uiz_varids[im][c])); - GRT_SAFE_ASPRINTF(&s_title, "r%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); - NC_CHECK(nc_inq_varid(in_ncid, s_title, &in_uir_varids[im][c])); + if(Ctrl->Q.active){ + // -Q：深度来自文件，禁止 -Dr（getopt 已拦一道，此处再保险） + if(Ctrl->D.r_active){ + GRTRaiseError("Do not set -Dr with -Q; receiver depths come from the points file."); } - GRT_SAFE_FREE_PTR(s_title); + } else { + // 网格接收：单台站深度库禁止 -Dr，多台站深度库必须 -Dr + if(!multi_rcv && Ctrl->D.r_active){ + GRTRaiseError("Library has a single receiver depth; do not set -Dr."); + } + if(multi_rcv && !Ctrl->D.r_active){ + GRTRaiseError("Library has multiple receiver depths; -Dr is required."); + } + } + + // 点源：多震源深度必须 -Ds（有限断层震源深度来自断层几何，不走此项） + if(Ctrl->isPointSource && multi_src && !Ctrl->D.s_active){ + GRTRaiseError("Library has multiple source depths; -Ds is required for point source."); + } +} + + +/** + * 构建接收点列表 + * + * -Q：任意点（各点自有深度，is_grid=false） + * 否则：-X/-Y 或延用库水平网格，统一深度（-Dr 或库 deprcvs[0]），is_grid=true + */ +static GRT_RECV_POINTS *build_syn_recv(const GRT_MODULE_CTRL *Ctrl, const STGRNLIB *lib) +{ + if(Ctrl->Q.active){ + return grt_recv_points_from_file(Ctrl->Q.s_path); } - // 定义合成结果 varid - for(int c=0; cisnewNEgrid ? Ctrl->X.nnorth : lib->nnorth; + size_t neast = Ctrl->isnewNEgrid ? Ctrl->Y.neast : lib->neast; + const real_t *norths = Ctrl->isnewNEgrid ? Ctrl->X.norths : lib->norths; + const real_t *easts = Ctrl->isnewNEgrid ? Ctrl->Y.easts : lib->easts; + real_t deprcv = Ctrl->D.r_active ? Ctrl->D.deprcv : lib->deprcvs[0]; + return grt_recv_points_from_grid(nnorth, norths, neast, easts, deprcv); +} + + +/** + * 定义位移分量及可选偏导变量 + * + * @param[in] ncid 输出 nc + * @param[in] ndims 维数（grid=2 / points=1） + * @param[in] dimids 维 id + * @param[in] chs 分量名字符（ZRT 或 ZNE） + * @param[in] calc_upar 是否定义偏导变量 + * @param[out] syn_varids + * @param[out] syn_upar_varids + */ +static void def_syn_channel_vars( + int ncid, int ndims, const int *dimids, const char *chs, bool calc_upar, + int syn_varids[GRT_CHANNEL_NUM], + int syn_upar_varids[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]) +{ + for(int c = 0; c < GRT_CHANNEL_NUM; ++c){ char *s_title = NULL; GRT_SAFE_ASPRINTF(&s_title, "%c", toupper(chs[c])); - NC_CHECK(nc_def_var(out_ncid, s_title, NC_REAL, ndims, out_dimids, &out_syn_varids[c])); - // 位移偏导 - if(Ctrl->e.active){ - for(int c2=0; c2depths[0]; + real_t rcv_va = 0.0, rcv_vb = 0.0, rcv_rho = 0.0; + grt_modarr_medium_at_depth( + lib->nlayer, lib->modarr, + deprcv, &rcv_va, &rcv_vb, &rcv_rho); + + NC_CHECK(NC_FUNC_REAL(nc_put_att)(ncid, NC_GLOBAL, "deprcv", NC_REAL, 1, &deprcv)); + NC_CHECK(NC_FUNC_REAL(nc_put_att)(ncid, NC_GLOBAL, "rcv_va", NC_REAL, 1, &rcv_va)); + NC_CHECK(NC_FUNC_REAL(nc_put_att)(ncid, NC_GLOBAL, "rcv_vb", NC_REAL, 1, &rcv_vb)); + NC_CHECK(NC_FUNC_REAL(nc_put_att)(ncid, NC_GLOBAL, "rcv_rho", NC_REAL, 1, &rcv_rho)); + + NC_CHECK(nc_def_dim(ncid, "north", recv->nnorth, &north_dimid)); + NC_CHECK(nc_def_dim(ncid, "east", recv->neast, &east_dimid)); + NC_CHECK(nc_def_var(ncid, "north", NC_REAL, 1, &north_dimid, &north_varid)); + NC_CHECK(nc_def_var(ncid, "east", NC_REAL, 1, &east_dimid, &east_varid)); + dimids[0] = north_dimid; + dimids[1] = east_dimid; + + def_syn_channel_vars(ncid, ndims, dimids, chs, calc_upar, syn_varids, syn_upar_varids); + NC_CHECK(nc_enddef(ncid)); + + // 从展开点列还原轴坐标（ipt = ieast + inorth*neast） + real_t *north_axis = (real_t *)calloc(recv->nnorth, sizeof(real_t)); + real_t *east_axis = (real_t *)calloc(recv->neast, sizeof(real_t)); + for(size_t inorth = 0; inorth < recv->nnorth; ++inorth){ + north_axis[inorth] = recv->norths[inorth * recv->neast]; + } + for(size_t ieast = 0; ieast < recv->neast; ++ieast){ + east_axis[ieast] = recv->easts[ieast]; + } + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, north_varid, north_axis)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, east_varid, east_axis)); + GRT_SAFE_FREE_PTR(north_axis); + GRT_SAFE_FREE_PTR(east_axis); +} - // 写入坐标变量 - NC_CHECK(NC_FUNC_REAL(nc_put_var) (out_ncid, out_x_varid, norths)); - NC_CHECK(NC_FUNC_REAL(nc_put_var) (out_ncid, out_y_varid, easts)); - // 总震中距数 - nr0 = nnorth0 * neast0; +/** + * 任意点布局：写逐点坐标，以及各点所在层的接收介质 + * + * 接收介质仅供后续应力使用，与点源/有限断层无关 + */ +static void put_syn_points_meta( + int ncid, const STGRNLIB *lib, const GRT_RECV_POINTS *recv, + int *syn_varids, int syn_upar_varids[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], + const char *chs, bool calc_upar) +{ + const int ndims = 1; + int point_dimid; + int north_varid, east_varid, depth_varid; + int rcv_va_varid, rcv_vb_varid, rcv_rho_varid; + int dimids[ndims]; + + NC_CHECK(nc_def_dim(ncid, "point", recv->npts, &point_dimid)); + dimids[0] = point_dimid; + + NC_CHECK(nc_def_var(ncid, "north", NC_REAL, ndims, dimids, &north_varid)); + NC_CHECK(nc_def_var(ncid, "east", NC_REAL, ndims, dimids, &east_varid)); + NC_CHECK(nc_def_var(ncid, "depth", NC_REAL, ndims, dimids, &depth_varid)); + NC_CHECK(nc_def_var(ncid, "rcv_va", NC_REAL, ndims, dimids, &rcv_va_varid)); + NC_CHECK(nc_def_var(ncid, "rcv_vb", NC_REAL, ndims, dimids, &rcv_vb_varid)); + NC_CHECK(nc_def_var(ncid, "rcv_rho", NC_REAL, ndims, dimids, &rcv_rho_varid)); + + def_syn_channel_vars(ncid, ndims, dimids, chs, calc_upar, syn_varids, syn_upar_varids); + NC_CHECK(nc_enddef(ncid)); + + real_t *rcv_va = (real_t *)calloc(recv->npts, sizeof(real_t)); + real_t *rcv_vb = (real_t *)calloc(recv->npts, sizeof(real_t)); + real_t *rcv_rho = (real_t *)calloc(recv->npts, sizeof(real_t)); + for(size_t i = 0; i < recv->npts; ++i){ + grt_modarr_medium_at_depth( + lib->nlayer, lib->modarr, + recv->depths[i], &rcv_va[i], &rcv_vb[i], &rcv_rho[i]); + } - // 先将所有格林函数及其偏导读入内存， - // 否则连续使用 nc_grt_var1 式读入效率太慢 - realChnlGrid *grn = (realChnlGrid *) calloc(nr0, sizeof(*grn)); - realChnlGrid *grn_uiz = (Ctrl->e.active)? (realChnlGrid *) calloc(nr0, sizeof(*grn_uiz)) : NULL; - realChnlGrid *grn_uir = (Ctrl->e.active)? (realChnlGrid *) calloc(nr0, sizeof(*grn_uir)) : NULL; - { - real_t *u = (real_t *)calloc(nr0, sizeof(real_t)); - GRT_LOOP_ChnlGrid(im, c){ - int modr = GRT_SRC_M_ORDERS[im]; - if(modr==0 && GRT_ZRT_CODES[c]=='T') continue; + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, north_varid, recv->norths)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, east_varid, recv->easts)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, depth_varid, recv->depths)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, rcv_va_varid, rcv_va)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, rcv_vb_varid, rcv_vb)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, rcv_rho_varid, rcv_rho)); - NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_u_varids[im][c], u)); - for(size_t ir = 0; ir < nr0; ++ir){ - grn[ir][im][c] = u[ir]; - } + GRT_SAFE_FREE_PTR(rcv_va); + GRT_SAFE_FREE_PTR(rcv_vb); + GRT_SAFE_FREE_PTR(rcv_rho); +} - if(Ctrl->e.active){ - NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_uiz_varids[im][c], u)); - for(size_t ir = 0; ir < nr0; ++ir){ - grn_uiz[ir][im][c] = u[ir]; - } - NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_uir_varids[im][c], u)); - for(size_t ir = 0; ir < nr0; ++ir){ - grn_uir[ir][im][c] = u[ir]; - } - } +/** + * 写入 syn / syn_upar 场量 + * + * @param[in] ncid + * @param[in] npts + * @param[in] calc_upar + * @param[in] syn + * @param[in] syn_upar + * @param[in] syn_varids + * @param[in] syn_upar_varids + */ +static void put_syn_fields( + int ncid, size_t npts, bool calc_upar, + const real_t (*syn)[GRT_CHANNEL_NUM], + const real_t (*syn_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM], + const int syn_varids[GRT_CHANNEL_NUM], + const int syn_upar_varids[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]) +{ + real_t *tmp = (real_t *)calloc(npts, sizeof(real_t)); + for(int c = 0; c < GRT_CHANNEL_NUM; ++c){ + for(size_t ir = 0; ir < npts; ++ir){ + tmp[ir] = syn[ir][c]; } - GRT_SAFE_FREE_PTR(u); - } - - // 新接收点网格总点数 - nr = nnorth * neast; - - // 最终计算的结果 - real_t (*syn)[GRT_CHANNEL_NUM] = (real_t (*)[GRT_CHANNEL_NUM])calloc(nr, sizeof(real_t)*GRT_CHANNEL_NUM); - real_t (*syn_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM] = (real_t (*)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM])calloc(nr, sizeof(real_t)*GRT_CHANNEL_NUM*GRT_CHANNEL_NUM); - - static_syn_from_gf( - nnorth0, norths0, neast0, easts0, - nnorth, norths, neast, easts, - grn, grn_uiz, grn_uir, - Ctrl->computeType, Ctrl->S.M0, src_va/src_vb, Ctrl->mchn, - rot2ZNE, Ctrl->e.active, - syn, syn_upar - ); - - // 写入 nc 文件 - real_t *tmpdata = (real_t *)calloc(nr, sizeof(real_t)); - for(int c=0; ce.active){ - for(int c2=0; c2N.active; + bool calc_upar = Ctrl->e.active; + const char *chs = rot2ZNE ? GRT_ZNE_CODES : GRT_ZRT_CODES; + const char *layout_str = recv->is_grid ? GRT_RECV_LAYOUT_GRID : GRT_RECV_LAYOUT_POINTS; + + int ncid; + int syn_varids[GRT_CHANNEL_NUM]; + int syn_upar_varids[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM]; + + NC_CHECK(nc_create(path, NC_CLOBBER, &ncid)); + NC_CHECK(nc_put_att_text(ncid, NC_GLOBAL, "layout", strlen(layout_str), layout_str)); + + // 点源才记录单一 depsrc；有限断层各子源深度不同，不写此属性 + if(Ctrl->isPointSource){ + NC_CHECK(NC_FUNC_REAL(nc_put_att)(ncid, NC_GLOBAL, "depsrc", NC_REAL, 1, &depsrc)); + } - // 关闭文件 - NC_CHECK(nc_close(in_ncid)); - NC_CHECK(nc_close(out_ncid)); + { + int int_calc_upar = calc_upar ? 1 : 0; + int int_rot2ZNE = rot2ZNE ? 1 : 0; + NC_CHECK(nc_put_att_int(ncid, NC_GLOBAL, "calc_upar", NC_INT, 1, &int_calc_upar)); + NC_CHECK(nc_put_att_int(ncid, NC_GLOBAL, "rot2ZNE", NC_INT, 1, &int_rot2ZNE)); + } + NC_CHECK(nc_put_att_text( + ncid, NC_GLOBAL, "computeType", + strlen(Ctrl->s_computeType), Ctrl->s_computeType)); + + // 坐标与接收介质：仅随布局变化，与点源/有限断层无关 + if(recv->is_grid){ + put_syn_grid_meta(ncid, lib, recv, syn_varids, syn_upar_varids, chs, calc_upar); + } else { + put_syn_points_meta(ncid, lib, recv, syn_varids, syn_upar_varids, chs, calc_upar); + } + + put_syn_fields(ncid, recv->npts, calc_upar, syn, syn_upar, syn_varids, syn_upar_varids); + NC_CHECK(nc_close(ncid)); +} - if(! Ctrl->s.active) { - GRTRaiseInfo("Synthetic static displacements of %s source saved in \"%s\".", srcTypeFullName[Ctrl->computeType], Ctrl->O.s_outgrid); + +/** 子模块主函数 */ +int static_syn_main(int argc, char **argv){ + GRT_MODULE_CTRL *Ctrl = calloc(1, sizeof(*Ctrl)); + getopt_from_command(Ctrl, argc, argv); + + STGRNLIB *lib = grt_stgrnlib_load_nc(Ctrl->G.s_ingrid); + check_syn_depth_options(Ctrl, lib); + if(Ctrl->e.active && !lib->calc_upar){ + GRTRaiseError("Input grid didn't have displacement derivatives, you can't set -e."); } - // 释放内存 - GRT_SAFE_FREE_PTR(norths0); - GRT_SAFE_FREE_PTR(easts0); + // 接收点：网格共面 或 -Q 逐点（深度范围在 PS/FF 内再校验） + GRT_RECV_POINTS *recv = build_syn_recv(Ctrl, lib); + + real_t (*syn)[GRT_CHANNEL_NUM] = (real_t (*)[GRT_CHANNEL_NUM])calloc( + recv->npts, sizeof(real_t) * GRT_CHANNEL_NUM); + real_t (*syn_upar)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM] = + (real_t (*)[GRT_CHANNEL_NUM][GRT_CHANNEL_NUM])calloc( + recv->npts, sizeof(real_t) * GRT_CHANNEL_NUM * GRT_CHANNEL_NUM); + + // depsrc 仅点源需要；有限断层由各子断层几何提供 + real_t depsrc = 0.0; + if(Ctrl->isPointSource){ + depsrc = Ctrl->D.s_active ? Ctrl->D.depsrc : lib->depsrcs[0]; + static_syn_from_gf_PS( + lib, depsrc, + recv->npts, recv->norths, recv->easts, recv->depths, + recv->is_grid, + Ctrl->computeType, Ctrl->S.M0, Ctrl->S.mult_src_mu, Ctrl->mchn, + Ctrl->N.active, Ctrl->e.active, + syn, syn_upar); + } else { + static_syn_from_gf_FF( + lib, + Ctrl->C.nfault, Ctrl->C.faults, + Ctrl->C.dL, Ctrl->C.dW, + recv->npts, recv->norths, recv->easts, recv->depths, + recv->is_grid, Ctrl->e.active, + syn, syn_upar); + } - GRT_SAFE_FREE_PTR(grn); - GRT_SAFE_FREE_PTR(grn_uiz); - GRT_SAFE_FREE_PTR(grn_uir); - GRT_SAFE_FREE_PTR(syn); - GRT_SAFE_FREE_PTR(syn_upar); + // 写出 nc（接收介质只在此处按布局查询，不参与合成） + save_syn_nc(Ctrl->O.s_outgrid, Ctrl, lib, recv, depsrc, syn, syn_upar); + if(!Ctrl->s.active){ + GRTRaiseInfo( + "Synthetic static displacements of %s source saved in \"%s\".", + srcTypeFullName[Ctrl->computeType], Ctrl->O.s_outgrid); + } + GRT_SAFE_FREE_PTR(syn); + GRT_SAFE_FREE_PTR(syn_upar); + grt_recv_points_free(recv); + grt_stgrnlib_free(lib); free_Ctrl(Ctrl); return EXIT_SUCCESS; } \ No newline at end of file diff --git a/pygrt/C_extension/src/static/recv_points.c b/pygrt/C_extension/src/static/recv_points.c new file mode 100644 index 00000000..45fa33a7 --- /dev/null +++ b/pygrt/C_extension/src/static/recv_points.c @@ -0,0 +1,149 @@ +/** + * @file recv_points.c + * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) + * @date 2026-08 + * + * 静态 syn / 后处理用的接收点列表 + * + */ + +#include +#include +#include + +#include "grt/static/recv_points.h" +#include "grt/common/checkerror.h" +#include "grt/common/util.h" +#include "grt/common/mynetcdf.h" + +GRT_RECV_POINTS *grt_recv_points_from_grid( + size_t nnorth, const real_t *norths, + size_t neast, const real_t *easts, + real_t depth) +{ + if(nnorth == 0 || neast == 0 || norths == NULL || easts == NULL){ + GRTRaiseError("empty receiver grid."); + } + if(depth < 0.0){ + GRTRaiseError("Negative receiver depth is not supported."); + } + + GRT_RECV_POINTS *pts = (GRT_RECV_POINTS *)calloc(1, sizeof(GRT_RECV_POINTS)); + pts->is_grid = true; + pts->nnorth = nnorth; + pts->neast = neast; + pts->npts = nnorth * neast; + pts->norths = (real_t *)calloc(pts->npts, sizeof(real_t)); + pts->easts = (real_t *)calloc(pts->npts, sizeof(real_t)); + pts->depths = (real_t *)calloc(pts->npts, sizeof(real_t)); + + for(size_t inorth = 0; inorth < nnorth; ++inorth){ + for(size_t ieast = 0; ieast < neast; ++ieast){ + size_t ipt = ieast + inorth * neast; + pts->norths[ipt] = norths[inorth]; + pts->easts[ipt] = easts[ieast]; + pts->depths[ipt] = depth; + } + } + return pts; +} + + +GRT_RECV_POINTS *grt_recv_points_from_file(const char *path) +{ + GRTCheckFileExist(path); + + FILE *fp = fopen(path, "r"); + if(fp == NULL){ + GRTRaiseError("Failed to open receiver points file \"%s\".", path); + } + + // 先统计有效行数 + size_t npts = 0; + char *line = NULL; + size_t nlen = 0; + while(grt_getline(&line, &nlen, fp) != -1){ + grt_trim_whitespace(line); + if(grt_is_comment_or_empty(line)) continue; + npts++; + } + if(npts == 0){ + GRT_SAFE_FREE_PTR(line); + fclose(fp); + GRTRaiseError("No receiver points found in \"%s\".", path); + } + + GRT_RECV_POINTS *pts = (GRT_RECV_POINTS *)calloc(1, sizeof(GRT_RECV_POINTS)); + pts->is_grid = false; + pts->nnorth = 0; + pts->neast = 0; + pts->npts = npts; + pts->norths = (real_t *)calloc(npts, sizeof(real_t)); + pts->easts = (real_t *)calloc(npts, sizeof(real_t)); + pts->depths = (real_t *)calloc(npts, sizeof(real_t)); + + rewind(fp); + size_t ipt = 0; + size_t lineno = 0; + while(grt_getline(&line, &nlen, fp) != -1){ + lineno++; + grt_trim_whitespace(line); + if(grt_is_comment_or_empty(line)) continue; + + real_t n, e, d; + if(3 != sscanf(line, "%lf %lf %lf", &n, &e, &d)){ + GRT_SAFE_FREE_PTR(line); + grt_recv_points_free(pts); + fclose(fp); + GRTRaiseError( + "Invalid receiver point at line %zu in \"%s\" (expect: north east depth).", + lineno, path); + } + if(d < 0.0){ + GRT_SAFE_FREE_PTR(line); + grt_recv_points_free(pts); + fclose(fp); + GRTRaiseError("Negative receiver depth at line %zu in \"%s\".", lineno, path); + } + pts->norths[ipt] = n; + pts->easts[ipt] = e; + pts->depths[ipt] = d; + ipt++; + } + + GRT_SAFE_FREE_PTR(line); + fclose(fp); + return pts; +} + + +void grt_recv_points_free(GRT_RECV_POINTS *pts) +{ + if(pts == NULL) return; + GRT_SAFE_FREE_PTR(pts->norths); + GRT_SAFE_FREE_PTR(pts->easts); + GRT_SAFE_FREE_PTR(pts->depths); + GRT_SAFE_FREE_PTR(pts); +} + + +bool grt_recv_nc_is_points(int ncid) +{ + size_t len = 0; + int status = nc_inq_attlen(ncid, NC_GLOBAL, "layout", &len); + if(status == NC_NOERR && len > 0){ + char *layout = (char *)calloc(len + 1, 1); + if(nc_get_att_text(ncid, NC_GLOBAL, "layout", layout) == NC_NOERR){ + bool is_pts = (strcmp(layout, GRT_RECV_LAYOUT_POINTS) == 0); + GRT_SAFE_FREE_PTR(layout); + return is_pts; + } + GRT_SAFE_FREE_PTR(layout); + } + + int dimid; + if(nc_inq_dimid(ncid, "point", &dimid) == NC_NOERR){ + return true; + } + return false; +} diff --git a/pygrt/C_extension/src/static/stgrnlib.c b/pygrt/C_extension/src/static/stgrnlib.c new file mode 100644 index 00000000..043d8088 --- /dev/null +++ b/pygrt/C_extension/src/static/stgrnlib.c @@ -0,0 +1,510 @@ +/** + * @file stgrnlib.c + * @author Zhu Dengda (zhudengda@mail.iggcas.ac.cn) + * @date 2026-08 + * + * 静态格林函数库 STGRNLIB：内存管理与四维 nc 读写 + * + */ + +#include +#include + +#include "grt/static/stgrnlib.h" +#include "grt/common/mynetcdf.h" +#include "grt/common/checkerror.h" +#include "grt/common/search.h" + +/** 检查数组严格升序 */ +static void require_strictly_ascending(const real_t *a, size_t n, const char *name) +{ + for(size_t i = 1; i < n; ++i){ + if(!(a[i] > a[i - 1])){ + GRTRaiseError("%s must be strictly ascending.", name); + } + } +} + +/** + * 由 norths/easts 填充网格序 rs，以及升序 sort_rs / sort_rs_idx / isUniform / dr + * + * isUniform 判定与 syn 一致：nr>2 且网格序已是等距升序 + */ +static void fill_rs_meta(STGRNLIB *lib) +{ + lib->nr = lib->nnorth * lib->neast; + lib->rs = (real_t *)malloc(lib->nr * sizeof(real_t)); + lib->sort_rs = (real_t *)malloc(lib->nr * sizeof(real_t)); + lib->sort_rs_idx = (size_t *)malloc(lib->nr * sizeof(size_t)); + + for(size_t inorth = 0; inorth < lib->nnorth; ++inorth){ + for(size_t ieast = 0; ieast < lib->neast; ++ieast){ + size_t ipt = ieast + inorth * lib->neast; + lib->rs[ipt] = hypot(lib->norths[inorth], lib->easts[ieast]); + lib->sort_rs_idx[ipt] = ipt; + } + } + memcpy(lib->sort_rs, lib->rs, lib->nr * sizeof(real_t)); + + // 还未排序前，先判断是否是一个等距升序数组，便于加快后续查找 + lib->isUniform = (lib->nr > 2); + lib->dr = (lib->nr > 1) ? (lib->rs[1] - lib->rs[0]) : 0.0; + for(size_t ir = 1; ir + 1 < lib->nr; ++ir){ + if(fabs(2.0 * lib->rs[ir] - (lib->rs[ir - 1] + lib->rs[ir + 1])) > 1e-3 + || lib->rs[ir - 1] >= lib->rs[ir] + || lib->rs[ir] >= lib->rs[ir + 1]){ + lib->isUniform = false; + break; + } + } + + if(!lib->isUniform){ + if(lib->nr > 1 + && grt_argsort(lib->rs, lib->nr, sizeof(*lib->rs), grt_compare_real_t, lib->sort_rs_idx) != 0){ + GRTRaiseError("Unable to sort epicentral distances."); + } + for(size_t i = 0; i < lib->nr; ++i){ + lib->sort_rs[i] = lib->rs[lib->sort_rs_idx[i]]; + } + } +} + +/** 按已填好的维度申请 u/uiz/uir */ +static void allocate_u(STGRNLIB *lib) +{ + size_t nr = lib->nr; + lib->u = (realChnlGrid ***)calloc(lib->ndepsrc, sizeof(*lib->u)); + lib->uiz = lib->calc_upar ? (realChnlGrid ***)calloc(lib->ndepsrc, sizeof(*lib->uiz)) : NULL; + lib->uir = lib->calc_upar ? (realChnlGrid ***)calloc(lib->ndepsrc, sizeof(*lib->uir)) : NULL; + + for(size_t is = 0; is < lib->ndepsrc; ++is){ + lib->u[is] = (realChnlGrid **)calloc(lib->ndeprcv, sizeof(*lib->u[is])); + if(lib->calc_upar){ + lib->uiz[is] = (realChnlGrid **)calloc(lib->ndeprcv, sizeof(*lib->uiz[is])); + lib->uir[is] = (realChnlGrid **)calloc(lib->ndeprcv, sizeof(*lib->uir[is])); + } + for(size_t ir = 0; ir < lib->ndeprcv; ++ir){ + lib->u[is][ir] = (realChnlGrid *)calloc(nr, sizeof(realChnlGrid)); + if(lib->calc_upar){ + lib->uiz[is][ir] = (realChnlGrid *)calloc(nr, sizeof(realChnlGrid)); + lib->uir[is][ir] = (realChnlGrid *)calloc(nr, sizeof(realChnlGrid)); + } + } + } +} + +/** 仅释放 u/uiz/uir */ +static void free_u(STGRNLIB *lib) +{ + if(lib == NULL) return; + if(lib->u != NULL){ + for(size_t is = 0; is < lib->ndepsrc; ++is){ + if(lib->u[is] != NULL){ + for(size_t ir = 0; ir < lib->ndeprcv; ++ir){ + GRT_SAFE_FREE_PTR(lib->u[is][ir]); + if(lib->uiz && lib->uiz[is]) GRT_SAFE_FREE_PTR(lib->uiz[is][ir]); + if(lib->uir && lib->uir[is]) GRT_SAFE_FREE_PTR(lib->uir[is][ir]); + } + } + GRT_SAFE_FREE_PTR(lib->u[is]); + if(lib->uiz) GRT_SAFE_FREE_PTR(lib->uiz[is]); + if(lib->uir) GRT_SAFE_FREE_PTR(lib->uir[is]); + } + } + GRT_SAFE_FREE_PTR(lib->u); + GRT_SAFE_FREE_PTR(lib->uiz); + GRT_SAFE_FREE_PTR(lib->uir); +} + +STGRNLIB *grt_stgrnlib_alloc( + size_t ndepsrc, const real_t *depsrcs, + size_t ndeprcv, const real_t *deprcvs, + size_t nnorth, const real_t *norths, + size_t neast, const real_t *easts, + bool calc_upar) +{ + if(ndepsrc == 0 || ndeprcv == 0 || nnorth == 0 || neast == 0){ + GRTRaiseError("dimensions must be positive."); + } + if(depsrcs == NULL || deprcvs == NULL || norths == NULL || easts == NULL){ + GRTRaiseError("coordinate arrays are NULL."); + } + require_strictly_ascending(depsrcs, ndepsrc, "depsrcs"); + require_strictly_ascending(deprcvs, ndeprcv, "deprcvs"); + require_strictly_ascending(norths, nnorth, "norths"); + require_strictly_ascending(easts, neast, "easts"); + + STGRNLIB *lib = (STGRNLIB *)calloc(1, sizeof(*lib)); + lib->ndepsrc = ndepsrc; + lib->ndeprcv = ndeprcv; + lib->nnorth = nnorth; + lib->neast = neast; + lib->calc_upar = calc_upar; + + lib->depsrcs = (real_t *)malloc(ndepsrc * sizeof(real_t)); + lib->deprcvs = (real_t *)malloc(ndeprcv * sizeof(real_t)); + lib->norths = (real_t *)malloc(nnorth * sizeof(real_t)); + lib->easts = (real_t *)malloc(neast * sizeof(real_t)); + memcpy(lib->depsrcs, depsrcs, ndepsrc * sizeof(real_t)); + memcpy(lib->deprcvs, deprcvs, ndeprcv * sizeof(real_t)); + memcpy(lib->norths, norths, nnorth * sizeof(real_t)); + memcpy(lib->easts, easts, neast * sizeof(real_t)); + + fill_rs_meta(lib); + + // 介质参数由调用方随后填入 + lib->src_va = (real_t *)calloc(ndepsrc, sizeof(real_t)); + lib->src_vb = (real_t *)calloc(ndepsrc, sizeof(real_t)); + lib->src_rho = (real_t *)calloc(ndepsrc, sizeof(real_t)); + lib->rcv_va = (real_t *)calloc(ndeprcv, sizeof(real_t)); + lib->rcv_vb = (real_t *)calloc(ndeprcv, sizeof(real_t)); + lib->rcv_rho = (real_t *)calloc(ndeprcv, sizeof(real_t)); + lib->nlayer = 0; + lib->modarr = NULL; + + allocate_u(lib); + return lib; +} + +void grt_stgrnlib_free(STGRNLIB *lib) +{ + if(lib == NULL) return; + free_u(lib); + GRT_SAFE_FREE_PTR(lib->depsrcs); + GRT_SAFE_FREE_PTR(lib->deprcvs); + GRT_SAFE_FREE_PTR(lib->norths); + GRT_SAFE_FREE_PTR(lib->easts); + GRT_SAFE_FREE_PTR(lib->rs); + GRT_SAFE_FREE_PTR(lib->sort_rs); + GRT_SAFE_FREE_PTR(lib->sort_rs_idx); + GRT_SAFE_FREE_PTR(lib->src_va); + GRT_SAFE_FREE_PTR(lib->src_vb); + GRT_SAFE_FREE_PTR(lib->src_rho); + GRT_SAFE_FREE_PTR(lib->rcv_va); + GRT_SAFE_FREE_PTR(lib->rcv_vb); + GRT_SAFE_FREE_PTR(lib->rcv_rho); + GRT_SAFE_FREE_PTR(lib->modarr); + free(lib); +} + +void grt_stgrnlib_set_modarr( + STGRNLIB *lib, size_t nlayer, const real_t (*modarr)[GRT_MODARR_NCOL]) +{ + if(lib == NULL){ + GRTRaiseError("lib is NULL."); + } + if(nlayer == 0 || modarr == NULL){ + GRTRaiseError("nlayer and modarr must be non-empty."); + } + GRT_SAFE_FREE_PTR(lib->modarr); + lib->modarr = (real_t (*)[GRT_MODARR_NCOL])malloc( + sizeof(real_t) * GRT_MODARR_NCOL * nlayer); + if(lib->modarr == NULL){ + GRTRaiseError("Failed to allocate modarr."); + } + memcpy(lib->modarr, modarr, sizeof(real_t) * GRT_MODARR_NCOL * nlayer); + lib->nlayer = nlayer; +} + +real_t grt_stgrnlib_default_subfault_size(const STGRNLIB *lib) +{ + const real_t atol = 1e-8; + + if(lib == NULL){ + GRTRaiseError("lib is NULL."); + } + if(lib->sort_rs == NULL){ + GRTRaiseError("STGRNLIB distance metadata is missing."); + } + + real_t dr = -1.0; + for(size_t i = 0; i + 1 < lib->nr; ++i){ + real_t d = lib->sort_rs[i + 1] - lib->sort_rs[i]; + if(d > atol && (dr < 0.0 || d < dr)) dr = d; + } + + real_t dz = -1.0; + for(size_t i = 0; i + 1 < lib->ndepsrc; ++i){ + real_t d = fabs(lib->depsrcs[i + 1] - lib->depsrcs[i]); + if(d > atol && (dz < 0.0 || d < dz)) dz = d; + } + + if(dr < 0.0 && dz < 0.0){ + GRTRaiseError( + "Cannot infer default subfault size: " + "need at least two distinct epicentral-distance or source-depth samples."); + } + if(dr < 0.0) return dz; + if(dz < 0.0) return dr; + return GRT_MIN(dr, dz); +} + +/** 读入一层 (is, ir) 的通道；ncid 已打开，变量为 4D */ +static void read_nc_channels_slice(STGRNLIB *lib, size_t is, size_t ir, int ncid) +{ + size_t nr = lib->nr; + real_t *buf = (real_t *)calloc(nr, sizeof(real_t)); + size_t start[4] = {is, ir, 0, 0}; + size_t count[4] = {1, 1, lib->nnorth, lib->neast}; + + GRT_LOOP_ChnlGrid(im, c){ + int modr = GRT_SRC_M_ORDERS[im]; + if(modr == 0 && GRT_ZRT_CODES[c] == 'T') continue; + + char *s_title = NULL; + int varid; + + GRT_SAFE_ASPRINTF(&s_title, "%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); + NC_CHECK(nc_inq_varid(ncid, s_title, &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_vara)(ncid, varid, start, count, buf)); + for(size_t ipt = 0; ipt < nr; ++ipt){ + lib->u[is][ir][ipt][im][c] = buf[ipt]; + } + GRT_SAFE_FREE_PTR(s_title); + + if(lib->calc_upar){ + GRT_SAFE_ASPRINTF(&s_title, "z%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); + NC_CHECK(nc_inq_varid(ncid, s_title, &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_vara)(ncid, varid, start, count, buf)); + for(size_t ipt = 0; ipt < nr; ++ipt){ + lib->uiz[is][ir][ipt][im][c] = buf[ipt]; + } + GRT_SAFE_FREE_PTR(s_title); + + GRT_SAFE_ASPRINTF(&s_title, "r%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); + NC_CHECK(nc_inq_varid(ncid, s_title, &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_vara)(ncid, varid, start, count, buf)); + for(size_t ipt = 0; ipt < nr; ++ipt){ + lib->uir[is][ir][ipt][im][c] = buf[ipt]; + } + GRT_SAFE_FREE_PTR(s_title); + } + } + + GRT_SAFE_FREE_PTR(buf); +} + +STGRNLIB *grt_stgrnlib_load_nc(const char *path) +{ + GRTCheckFileExist(path); + + int ncid; + NC_CHECK(nc_open(path, NC_NOWRITE, &ncid)); + + int depsrc_dimid, deprcv_dimid, north_dimid, east_dimid; + size_t ndepsrc, ndeprcv, nnorth, neast; + NC_CHECK(nc_inq_dimid(ncid, "depsrc", &depsrc_dimid)); + NC_CHECK(nc_inq_dimlen(ncid, depsrc_dimid, &ndepsrc)); + NC_CHECK(nc_inq_dimid(ncid, "deprcv", &deprcv_dimid)); + NC_CHECK(nc_inq_dimlen(ncid, deprcv_dimid, &ndeprcv)); + NC_CHECK(nc_inq_dimid(ncid, "north", &north_dimid)); + NC_CHECK(nc_inq_dimlen(ncid, north_dimid, &nnorth)); + NC_CHECK(nc_inq_dimid(ncid, "east", &east_dimid)); + NC_CHECK(nc_inq_dimlen(ncid, east_dimid, &neast)); + + if(ndepsrc == 0 || ndeprcv == 0 || nnorth == 0 || neast == 0){ + NC_CHECK(nc_close(ncid)); + GRTRaiseError("Invalid STGRNLIB nc \"%s\": empty dimension.", path); + } + + int int_calc_upar = 0; + NC_CHECK(nc_get_att_int(ncid, NC_GLOBAL, "calc_upar", &int_calc_upar)); + + // 坐标轴读入后经 grt_stgrnlib_alloc 建壳 + real_t *depsrcs = (real_t *)calloc(ndepsrc, sizeof(real_t)); + real_t *deprcvs = (real_t *)calloc(ndeprcv, sizeof(real_t)); + real_t *norths = (real_t *)calloc(nnorth, sizeof(real_t)); + real_t *easts = (real_t *)calloc(neast, sizeof(real_t)); + + int varid; + NC_CHECK(nc_inq_varid(ncid, "depsrc", &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var)(ncid, varid, depsrcs)); + NC_CHECK(nc_inq_varid(ncid, "deprcv", &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var)(ncid, varid, deprcvs)); + NC_CHECK(nc_inq_varid(ncid, "north", &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var)(ncid, varid, norths)); + NC_CHECK(nc_inq_varid(ncid, "east", &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var)(ncid, varid, easts)); + + STGRNLIB *lib = grt_stgrnlib_alloc( + ndepsrc, depsrcs, ndeprcv, deprcvs, nnorth, norths, neast, easts, (int_calc_upar != 0)); + GRT_SAFE_FREE_PTR(depsrcs); + GRT_SAFE_FREE_PTR(deprcvs); + GRT_SAFE_FREE_PTR(norths); + GRT_SAFE_FREE_PTR(easts); + + NC_CHECK(nc_inq_varid(ncid, "src_va", &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var)(ncid, varid, lib->src_va)); + NC_CHECK(nc_inq_varid(ncid, "src_vb", &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var)(ncid, varid, lib->src_vb)); + NC_CHECK(nc_inq_varid(ncid, "src_rho", &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var)(ncid, varid, lib->src_rho)); + NC_CHECK(nc_inq_varid(ncid, "rcv_va", &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var)(ncid, varid, lib->rcv_va)); + NC_CHECK(nc_inq_varid(ncid, "rcv_vb", &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var)(ncid, varid, lib->rcv_vb)); + NC_CHECK(nc_inq_varid(ncid, "rcv_rho", &varid)); + NC_CHECK(NC_FUNC_REAL(nc_get_var)(ncid, varid, lib->rcv_rho)); + + // 建库模型矩阵（必需） + { + int layer_dimid, param_dimid; + size_t nlayer = 0, nparam = 0; + if(nc_inq_dimid(ncid, "layer", &layer_dimid) != NC_NOERR || + nc_inq_dimid(ncid, "model_param", ¶m_dimid) != NC_NOERR){ + NC_CHECK(nc_close(ncid)); + grt_stgrnlib_free(lib); + GRTRaiseError( + "Invalid STGRNLIB nc \"%s\": missing model dimensions " + "(layer, model_param); rebuild with current greenfn.", + path); + } + NC_CHECK(nc_inq_dimlen(ncid, layer_dimid, &nlayer)); + NC_CHECK(nc_inq_dimlen(ncid, param_dimid, &nparam)); + if(nlayer == 0 || nparam != GRT_MODARR_NCOL){ + NC_CHECK(nc_close(ncid)); + grt_stgrnlib_free(lib); + GRTRaiseError( + "Invalid STGRNLIB nc \"%s\": layer=%zu, model_param=%zu " + "(expect layer>0, model_param=%d).", + path, nlayer, nparam, GRT_MODARR_NCOL); + } + NC_CHECK(nc_inq_varid(ncid, "model", &varid)); + real_t (*modarr)[GRT_MODARR_NCOL] = (real_t (*)[GRT_MODARR_NCOL])malloc( + sizeof(real_t) * GRT_MODARR_NCOL * nlayer); + NC_CHECK(NC_FUNC_REAL(nc_get_var)(ncid, varid, (real_t *)modarr)); + grt_stgrnlib_set_modarr(lib, nlayer, (const real_t (*)[GRT_MODARR_NCOL])modarr); + GRT_SAFE_FREE_PTR(modarr); + } + + for(size_t is = 0; is < ndepsrc; ++is){ + for(size_t ir = 0; ir < ndeprcv; ++ir){ + read_nc_channels_slice(lib, is, ir, ncid); + } + } + + NC_CHECK(nc_close(ncid)); + return lib; +} + +void grt_stgrnlib_save_nc(const STGRNLIB *lib, const char *path) +{ + if(lib == NULL || path == NULL){ + GRTRaiseError("lib/path is NULL."); + } + if(lib->ndepsrc == 0 || lib->ndeprcv == 0 || lib->nnorth == 0 || lib->neast == 0){ + GRTRaiseError("empty library."); + } + if(lib->nlayer == 0 || lib->modarr == NULL){ + GRTRaiseError("STGRNLIB has no model matrix; call grt_stgrnlib_set_modarr first."); + } + + int ncid; + int depsrc_dimid, deprcv_dimid, north_dimid, east_dimid; + int layer_dimid, param_dimid; + int dimids[4]; + int depsrc_varid, deprcv_varid, north_varid, east_varid; + int model_varid; + int src_va_varid, src_vb_varid, src_rho_varid; + int rcv_va_varid, rcv_vb_varid, rcv_rho_varid; + intChnlGrid u_varids; + intChnlGrid uiz_varids; + intChnlGrid uir_varids; + + NC_CHECK(nc_create(path, NC_CLOBBER, &ncid)); + + { + int tmp = lib->calc_upar ? 1 : 0; + NC_CHECK(nc_put_att_int(ncid, NC_GLOBAL, "calc_upar", NC_INT, 1, &tmp)); + } + + NC_CHECK(nc_def_dim(ncid, "depsrc", lib->ndepsrc, &depsrc_dimid)); + NC_CHECK(nc_def_dim(ncid, "deprcv", lib->ndeprcv, &deprcv_dimid)); + NC_CHECK(nc_def_dim(ncid, "north", lib->nnorth, &north_dimid)); + NC_CHECK(nc_def_dim(ncid, "east", lib->neast, &east_dimid)); + NC_CHECK(nc_def_dim(ncid, "layer", lib->nlayer, &layer_dimid)); + NC_CHECK(nc_def_dim(ncid, "model_param", GRT_MODARR_NCOL, ¶m_dimid)); + dimids[0] = depsrc_dimid; + dimids[1] = deprcv_dimid; + dimids[2] = north_dimid; + dimids[3] = east_dimid; + + { + int model_dimids[2] = {layer_dimid, param_dimid}; + NC_CHECK(nc_def_var(ncid, "model", NC_REAL, 2, model_dimids, &model_varid)); + } + + NC_CHECK(nc_def_var(ncid, "depsrc", NC_REAL, 1, &depsrc_dimid, &depsrc_varid)); + NC_CHECK(nc_def_var(ncid, "deprcv", NC_REAL, 1, &deprcv_dimid, &deprcv_varid)); + NC_CHECK(nc_def_var(ncid, "north", NC_REAL, 1, &north_dimid, &north_varid)); + NC_CHECK(nc_def_var(ncid, "east", NC_REAL, 1, &east_dimid, &east_varid)); + + NC_CHECK(nc_def_var(ncid, "src_va", NC_REAL, 1, &depsrc_dimid, &src_va_varid)); + NC_CHECK(nc_def_var(ncid, "src_vb", NC_REAL, 1, &depsrc_dimid, &src_vb_varid)); + NC_CHECK(nc_def_var(ncid, "src_rho", NC_REAL, 1, &depsrc_dimid, &src_rho_varid)); + NC_CHECK(nc_def_var(ncid, "rcv_va", NC_REAL, 1, &deprcv_dimid, &rcv_va_varid)); + NC_CHECK(nc_def_var(ncid, "rcv_vb", NC_REAL, 1, &deprcv_dimid, &rcv_vb_varid)); + NC_CHECK(nc_def_var(ncid, "rcv_rho", NC_REAL, 1, &deprcv_dimid, &rcv_rho_varid)); + + GRT_LOOP_ChnlGrid(im, c){ + int modr = GRT_SRC_M_ORDERS[im]; + char *s_title = NULL; + if(modr == 0 && GRT_ZRT_CODES[c] == 'T') continue; + + GRT_SAFE_ASPRINTF(&s_title, "%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); + NC_CHECK(nc_def_var(ncid, s_title, NC_REAL, 4, dimids, &u_varids[im][c])); + if(lib->calc_upar){ + GRT_SAFE_ASPRINTF(&s_title, "z%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); + NC_CHECK(nc_def_var(ncid, s_title, NC_REAL, 4, dimids, &uiz_varids[im][c])); + GRT_SAFE_ASPRINTF(&s_title, "r%s%c", GRT_SRC_M_NAME_ABBR[im], GRT_ZRT_CODES[c]); + NC_CHECK(nc_def_var(ncid, s_title, NC_REAL, 4, dimids, &uir_varids[im][c])); + } + GRT_SAFE_FREE_PTR(s_title); + } + + NC_CHECK(nc_enddef(ncid)); + + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, depsrc_varid, lib->depsrcs)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, deprcv_varid, lib->deprcvs)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, north_varid, lib->norths)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, east_varid, lib->easts)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, src_va_varid, lib->src_va)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, src_vb_varid, lib->src_vb)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, src_rho_varid, lib->src_rho)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, rcv_va_varid, lib->rcv_va)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, rcv_vb_varid, lib->rcv_vb)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, rcv_rho_varid, lib->rcv_rho)); + NC_CHECK(NC_FUNC_REAL(nc_put_var)(ncid, model_varid, (const real_t *)lib->modarr)); + + size_t nr = lib->nr; + real_t *tmpdata = (real_t *)calloc(nr, sizeof(real_t)); + for(size_t is = 0; is < lib->ndepsrc; ++is){ + for(size_t ir = 0; ir < lib->ndeprcv; ++ir){ + size_t start[4] = {is, ir, 0, 0}; + size_t count[4] = {1, 1, lib->nnorth, lib->neast}; + + GRT_LOOP_ChnlGrid(im, c){ + int modr = GRT_SRC_M_ORDERS[im]; + if(modr == 0 && GRT_ZRT_CODES[c] == 'T') continue; + + for(size_t ipt = 0; ipt < nr; ++ipt){ + tmpdata[ipt] = lib->u[is][ir][ipt][im][c]; + } + NC_CHECK(NC_FUNC_REAL(nc_put_vara)(ncid, u_varids[im][c], start, count, tmpdata)); + + if(lib->calc_upar){ + for(size_t ipt = 0; ipt < nr; ++ipt){ + tmpdata[ipt] = lib->uiz[is][ir][ipt][im][c]; + } + NC_CHECK(NC_FUNC_REAL(nc_put_vara)(ncid, uiz_varids[im][c], start, count, tmpdata)); + for(size_t ipt = 0; ipt < nr; ++ipt){ + tmpdata[ipt] = lib->uir[is][ir][ipt][im][c]; + } + NC_CHECK(NC_FUNC_REAL(nc_put_vara)(ncid, uir_varids[im][c], start, count, tmpdata)); + } + } + } + } + GRT_SAFE_FREE_PTR(tmpdata); + NC_CHECK(nc_close(ncid)); +} diff --git a/pygrt/pymod.py b/pygrt/pymod.py index d2ed159c..f0f5a496 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -10,6 +10,7 @@ from __future__ import annotations import os +import warnings from ctypes import c_size_t, cast, c_void_p from pathlib import Path from typing import Dict, Iterable, Optional, Sequence, Union @@ -24,6 +25,7 @@ PathLike = Union[str, os.PathLike] +DepthLike = Union[float, Sequence[float]] __all__ = ["PyModel1D"] @@ -44,6 +46,33 @@ def _normalize_distarr(distarr): raise ValueError("distarr must be a scalar or a 1-D sequence of floats.") +def _normalize_depths(depths: DepthLike, name: str) -> np.ndarray: + """ + 将震源/接收深度规范为一维 float64 数组 + + 接受标量或一维浮点序列；空数组与负深度会报错 + """ + if isinstance(depths, (str, bytes)): + raise TypeError(f"{name} must be a float or a 1-D sequence of floats, not a string.") + arr = np.asarray(depths, dtype=np.float64) + if arr.ndim == 0: + arr = np.ascontiguousarray([float(arr)], dtype=np.float64) + elif arr.ndim == 1: + arr = np.ascontiguousarray(arr, dtype=np.float64) + else: + raise ValueError(f"{name} must be a scalar or a 1-D sequence of floats.") + if arr.size == 0: + raise ValueError(f"{name} must not be empty.") + if np.any(arr < 0.0): + raise ValueError(f"{name} must be nonnegative.") + return arr + + +def _format_depth_list(depths: np.ndarray) -> str: + """将深度数组格式化为 CLI ``-Ds``/``-Dr`` 的逗号列表""" + return ",".join(format_float(float(z)) for z in depths) + + class PyModel1D: """ File-based 1D layered model for GRT calculations. @@ -350,8 +379,8 @@ def compute_grn( def compute_static_grn( self, *, - depsrc: float, - deprcv: float, + depsrc: DepthLike, + deprcv: DepthLike, norths: Optional[Sequence[float]] = None, easts: Optional[Sequence[float]] = None, distarr: Optional[Sequence[float]] = None, @@ -370,8 +399,14 @@ def compute_static_grn( Compute static Green's functions with the ``grt static greenfn`` command. Call :meth:`set_static_grn_path` first. Results are written to the - configured NetCDF file and currently overwrite any existing content. - All arguments must be passed by keyword. + configured NetCDF file (4D STGRNLIB layout + ``[depsrc][deprcv][north][east]``) and currently overwrite any existing + content. All arguments must be passed by keyword. + + ``depsrc`` / ``deprcv`` may be a scalar or a 1-D sequence: + + * Single depth pair: CLI ``-Ddepsrc/deprcv``. + * Multiple depths: CLI ``-Ds...`` / ``-Dr...`` (comma-separated list). Receiver locations can be specified in either of two ways: @@ -380,8 +415,10 @@ def compute_static_grn( 2. ``distarr``, a list of epicentral distances in km. This is equivalent to placing receivers along the east axis with north = 0. - :param depsrc: Source depth in km. - :param deprcv: Receiver depth in km. + :param depsrc: Source depth(s) in km. Multiple values must be + strictly ascending. + :param deprcv: Receiver depth(s) in km. Multiple values must be + strictly ascending. :param norths: Three values defining the north-coordinate option ``-Xstart/stop/step`` in km. :param easts: Three values defining the east-coordinate @@ -414,19 +451,30 @@ def compute_static_grn( displacement. Required later if strain, stress or rotation will be computed. :param stats: Whether to write integration statistics. + Only available for a single source/receiver depth; + ignored with a warning for multi-depth runs. :return: ``None``. Results are written to the configured NetCDF file. """ if self.static_grn_path is None: raise RuntimeError("Call set_static_grn_path() before compute_static_grn().") - if depsrc < 0 or deprcv < 0: - raise ValueError("Source and receiver depths must be nonnegative.") + + depsrcs = _normalize_depths(depsrc, "depsrc") + deprcvs = _normalize_depths(deprcv, "deprcv") + if depsrcs.size > 1 and not np.all(np.diff(depsrcs) > 0.0): + raise ValueError("depsrc must be strictly ascending when multiple values are given.") + if deprcvs.size > 1 and not np.all(np.diff(deprcvs) > 0.0): + raise ValueError("deprcv must be strictly ascending when multiple values are given.") + multi_depth = (depsrcs.size > 1) or (deprcvs.size > 1) + if distarr is not None: if norths is not None or easts is not None: raise ValueError("Use either distarr or norths/easts.") _, distances = _normalize_distarr(distarr) if distances.size == 0 or np.any(distances < 0.0): raise ValueError("distarr must contain nonnegative distances.") + if distances.size > 1 and not np.all(np.diff(distances) > 0.0): + raise ValueError("distarr must be strictly ascending.") command_grid = { "R": f"-R{','.join(format_float(value) for value in distances)}" } @@ -442,10 +490,14 @@ def compute_static_grn( "module": "static", "subcommand": "greenfn", "M": f"-M{self.modelpath}", - "D": f"-D{format_float(depsrc)}/{format_float(deprcv)}", - "O": f"-O{self.static_grn_path}", - "B": f"-B{self._boundary_option()}", } + if multi_depth: + command["Ds"] = f"-Ds{_format_depth_list(depsrcs)}" + command["Dr"] = f"-Dr{_format_depth_list(deprcvs)}" + else: + command["D"] = f"-D{format_float(float(depsrcs[0]))}/{format_float(float(deprcvs[0]))}" + command["O"] = f"-O{self.static_grn_path}" + command["B"] = f"-B{self._boundary_option()}" command.update(command_grid) # Build the -L option. @@ -472,7 +524,13 @@ def compute_static_grn( # Build the statistics and derivative options. if stats: - command["S"] = "-S" + if multi_depth: + warnings.warn( + "stats is ignored for multi-depth STGRNLIB computation.", + stacklevel=2, + ) + else: + command["S"] = "-S" if calc_upar: command["e"] = "-e" @@ -618,8 +676,8 @@ def compute_syn( def compute_static_syn( self, *, - scale: float, output_path: PathLike, + scale: Optional[float] = None, source: str = "EX", strike: Optional[float] = None, dip: Optional[float] = None, @@ -627,8 +685,13 @@ def compute_static_syn( force: Optional[Sequence[float]] = None, moment_tensor: Optional[Sequence[float]] = None, scale_with_mu: bool = False, + depsrc: Optional[float] = None, + deprcv: Optional[float] = None, norths: Optional[Sequence[float]] = None, easts: Optional[Sequence[float]] = None, + recv_points: Optional[PathLike] = None, + finite_fault: Optional[PathLike] = None, + subfault_size: Optional[Sequence[float]] = None, zne: bool = False, calc_upar: bool = False, return_result: bool = False, @@ -636,17 +699,25 @@ def compute_static_syn( r""" Synthesize static three-component displacement with ``grt static syn``. - Results are written to the NetCDF file ``output_path``. Source-type and - component conventions match :meth:`compute_syn`. Call + Results are written to the NetCDF file ``output_path``. Call :meth:`set_static_grn_path` and :meth:`compute_static_grn` first. All arguments must be passed by keyword. - By default the output grid inherits the north/east grid of the static - Green's function file. You may pass ``norths`` and ``easts`` to request a - new grid; each node then uses the nearest epicentral-distance Green's - function, which is an approximation that reuses an existing library. + Receivers default to the library north/east grid. Optionally redefine + them with ``norths``/``easts`` (uniform ``deprcv`` when the library has + multiple receiver depths), or with ``recv_points`` for an ASCII file of + arbitrary ``north east depth`` points (CLI ``-Q``). ``recv_points`` is + mutually exclusive with ``norths``/``easts`` and ``deprcv``. If the + library was built with ``distarr`` / ``-R``, the default grid is a 1-D + line (north = 0, east = R); set ``norths``/``easts`` or ``recv_points`` + to obtain a 2-D field. - Choose one source type with ``source``: + Point sources use ``scale`` and ``source``. Finite faults use + ``finite_fault`` (Coulomb-format file, CLI ``-C``) instead; that path + requires a multi-source-depth library, automatically writes ZNE, and + ignores point-source options. + + Choose one point-source type with ``source``: * ``EX`` - explosion. Only ``scale`` is required. * ``DC`` - double-couple / shear. Requires ``strike``, ``dip`` and ``rake``. @@ -655,15 +726,18 @@ def compute_static_syn( * ``MT`` - moment tensor. Requires ``moment_tensor=(Mxx, Mxy, Mxz, Myy, Myz, Mzz)``. - :param scale: Source scaling factor. For ``EX``, ``DC``, + :param output_path: Output NetCDF file path. + :param scale: Point-source scaling factor. For ``EX``, ``DC``, ``TS`` and ``MT``, this is the scalar seismic moment in dyne·cm. For ``SF``, the unit is dyne. If ``scale_with_mu`` is true, ``scale`` is treated as area × slip in cm³ and multiplied by the source-layer shear modulus :math:`\mu`. - :param output_path: Output NetCDF file path. - :param source: Source type. One of ``EX``, ``DC``, ``TS``, - ``SF`` and ``MT``. + Required for point sources; ignored for + ``finite_fault``. + :param source: Point-source type. One of ``EX``, ``DC``, + ``TS``, ``SF`` and ``MT``. Ignored when + ``finite_fault`` is set. :param strike: Fault strike in deg, in [0, 360]. North is 0°, clockwise positive. Required for ``DC`` and ``TS``. @@ -680,18 +754,37 @@ def compute_static_syn( Subscripts x/y/z denote north/east/down. :param scale_with_mu: If true, multiply ``scale`` by the source-layer shear modulus :math:`\mu` (CLI ``-Su``). + :param depsrc: Point-source depth in km (CLI ``-Ds``). Required + when the library has multiple source depths; + forbidden for ``finite_fault``. + :param deprcv: Receiver depth in km for grid receivers + (CLI ``-Dr``). Required when the library has + multiple receiver depths and ``recv_points`` + is not used. :param norths: Optional new north grid as three values ``(start, stop, step)`` in km. Must be set - together with ``easts``. + together with ``easts``. Mutually exclusive + with ``recv_points``. :param easts: Optional new east grid as three values ``(start, stop, step)`` in km. Must be set - together with ``norths``. + together with ``norths``. Mutually exclusive + with ``recv_points``. + :param recv_points: ASCII file of arbitrary receivers + (``north east depth`` in km; ``#`` comments). + Mutually exclusive with ``norths``/``easts`` + and ``deprcv``. + :param finite_fault: Coulomb-format finite-fault file (CLI ``-C``). + Mutually exclusive with point-source options. + :param subfault_size: Optional ``(dL, dW)`` in km for finite-fault + subdivision along strike / dip. If omitted, + the C code uses ``min(dr, dz)`` of the library. :param zne: If true, output ZNE instead of ZRT components. + Finite faults always write ZNE. :param calc_upar: If true, also synthesize spatial derivatives of - displacement. Derivative variable names are - prefixed with ``z``, ``r`` or ``t``. Set this - when strain, stress or rotation will be computed - later. + displacement. Derivative variable names use + prefixes ``z``/``r``/``t`` (ZRT) or + ``z``/``n``/``e`` (ZNE). Set this when strain, + stress or rotation will be computed later. :param return_result: If true, read the generated NetCDF file with :func:`pygrt.utils.read_static_nc`. @@ -702,25 +795,69 @@ def compute_static_syn( raise RuntimeError("Call set_static_grn_path() before compute_static_syn().") output = Path(output_path) output.parent.mkdir(parents=True, exist_ok=True) + + use_ff = finite_fault is not None + use_q = recv_points is not None + use_xy = norths is not None or easts is not None + if use_ff and ( + scale is not None + or force is not None + or moment_tensor is not None + or strike is not None + or dip is not None + or rake is not None + or scale_with_mu + or source.upper() != "EX" + ): + raise ValueError("finite_fault is mutually exclusive with point-source options.") + if use_q and use_xy: + raise ValueError("recv_points is mutually exclusive with norths/easts.") + if use_q and deprcv is not None: + raise ValueError("recv_points is mutually exclusive with deprcv.") + if use_xy and (norths is None or easts is None): + raise ValueError("norths and easts must be supplied together.") + if depsrc is not None and depsrc < 0.0: + raise ValueError("depsrc must be nonnegative.") + if deprcv is not None and deprcv < 0.0: + raise ValueError("deprcv must be nonnegative.") + if use_ff and depsrc is not None: + raise ValueError("depsrc is forbidden when finite_fault is set.") + if subfault_size is not None: + if not use_ff: + raise ValueError("subfault_size requires finite_fault.") + if len(subfault_size) != 2: + raise ValueError("subfault_size must be (dL, dW).") + command = { "module": "static", "subcommand": "syn", "G": f"-G{self.static_grn_path}", - "S": f"-S{'u' if scale_with_mu else ''}{format_float(scale)}", "O": f"-O{output}", } - command.update( - self._source_options(source, strike, dip, rake, force, moment_tensor) - ) - # Build the coordinate options. - if norths is not None or easts is not None: - if norths is None or easts is None: - raise ValueError("norths and easts must be supplied together.") + if use_ff: + c_opt = f"-C{Path(finite_fault)}" + if subfault_size is not None: + c_opt += f"+i{format_float(subfault_size[0])}/{format_float(subfault_size[1])}" + command["C"] = c_opt + else: + if scale is None: + raise ValueError("scale is required for point-source synthesis.") + command["S"] = f"-S{'u' if scale_with_mu else ''}{format_float(scale)}" + command.update( + self._source_options(source, strike, dip, rake, force, moment_tensor) + ) + if depsrc is not None: + command["Ds"] = f"-Ds{format_float(depsrc)}" + + if deprcv is not None: + command["Dr"] = f"-Dr{format_float(deprcv)}" + if use_q: + command["Q"] = f"-Q{Path(recv_points)}" + elif use_xy: command["X"] = f"-X{format_range(norths, 'norths')}" command["Y"] = f"-Y{format_range(easts, 'easts')}" - # Build the component and derivative options. if zne: command["N"] = "-N" if calc_upar: diff --git a/test/_compare_c_py/compare_func.py b/test/_compare_c_py/compare_func.py index e8489c02..76f1afa0 100644 --- a/test/_compare_c_py/compare_func.py +++ b/test/_compare_c_py/compare_func.py @@ -75,11 +75,11 @@ def compare_sac_dirs( def _nc_variable_map(path: PathLike) -> dict: - """读取 NetCDF 变量数据，跳过坐标轴 north/east""" + """读取 NetCDF 变量数据，跳过坐标轴 north/east/depsrc/deprcv/depth""" result = {} with netcdf_file(str(path), mmap=False) as dataset: for name, variable in dataset.variables.items(): - if name in {"north", "east"}: + if name in {"north", "east", "depsrc", "deprcv", "depth"}: continue result[name] = np.array(variable[:], copy=True) return result diff --git a/test/_compare_c_py/compare_stgrnlib.py b/test/_compare_c_py/compare_stgrnlib.py new file mode 100644 index 00000000..6d82bdab --- /dev/null +++ b/test/_compare_c_py/compare_stgrnlib.py @@ -0,0 +1,313 @@ +""" +Compare STGRNLIB nc results: +1) CLI vs Python for single/multi source and receiver depths +2) Multi-depth library slices vs corresponding single-depth results +3) Synthesis from multi-depth libraries: CLI vs Python, exact-depth vs + single-depth library, depth-interpolation linearity, and tensor postprocess +""" + +from pathlib import Path + +import numpy as np +import pygrt +from pygrt.cli import format_float, run_grt +from scipy.io import netcdf_file + +from compare_func import compare_nc_files + +# 与 shell 脚本保持一致 +NORTHS = (-2.0, 2.0, 1.0) +EASTS = (-2.0, 2.0, 1.0) +DEPSRCS = np.array([1.0, 2.0, 3.0]) +DEPRCVS = np.array([0.0, 0.5, 1.0]) + +MODNAME = str((Path(__file__).resolve().parent.parent / "milrow")) +CMPDIR = Path("stgrnlib_cmp") +ATOL = 1e-10 +RTOL = 1e-8 +# 平均相对误差阈值 +MEAN_RERR_MAX = 1e-6 +SCALE = 1e20 + +SYN_SKIP = { + "north", "east", "depth", "depsrc", "deprcv", + "rcv_va", "rcv_vb", "rcv_rho", "model", +} + + +def _coord_tag(z: float) -> str: + """与 shell 脚本一致：整数写 0/1，小数点换成 p（如 0p5）""" + zf = float(z) + if zf.is_integer(): + return str(int(zf)) + return str(zf).replace(".", "p") + + +def load_stgrnlib_nc(path: str) -> dict: + """读取 STGRNLIB nc 为字典（坐标 + 各通道数组）""" + out = {} + with netcdf_file(path, mmap=False) as f: + for k, v in f.variables.items(): + out[k] = np.array(v[:], dtype=float).copy() + return out + + +def compare_nc_data(a: dict, b: dict, label: str) -> float: + """逐通道比较，返回平均相对误差""" + for key in ("depsrc", "deprcv", "north", "east"): + if not np.allclose(a[key], b[key], rtol=RTOL, atol=ATOL): + raise ValueError(f"{label}: coordinate mismatch on '{key}'\n" + f" a={a[key]}\n b={b[key]}") + + med_keys = ("src_va", "src_vb", "src_rho", "rcv_va", "rcv_vb", "rcv_rho") + for key in med_keys: + if key in a and key in b: + if not np.allclose(a[key], b[key], rtol=RTOL, atol=ATOL): + raise ValueError(f"{label}: medium mismatch on '{key}'\n" + f" a={a[key]}\n b={b[key]}") + + skip = {"depsrc", "deprcv", "north", "east", "model", *med_keys} + keys = sorted(set(a.keys()) & set(b.keys()) - skip) + if not keys: + raise ValueError(f"{label}: no channel variables to compare") + + errors = [] + for k in keys: + va, vb = a[k], b[k] + if va.ndim != 4 or vb.ndim != 4: + continue + if va.shape != vb.shape: + raise ValueError(f"{label}: shape mismatch on '{k}': {va.shape} vs {vb.shape}") + if np.all(va == 0.0) and np.all(vb == 0.0): + continue + denom = np.mean(np.abs(vb)) + if denom == 0.0: + denom = np.mean(np.abs(va)) + rerr = np.sum(np.abs(va - vb)) / denom + print(f" {k}: {rerr:.6e}") + errors.append(rerr) + + mean_err = float(np.mean(errors)) if errors else 0.0 + print(f"{label}: mean relative error = {mean_err:.6e}") + if mean_err > MEAN_RERR_MAX: + raise ValueError(f"{label}: mean relative error too large ({mean_err})") + return mean_err + + +def extract_slice(lib: dict, isrc: int, ircv: int) -> dict: + """从多深度库取出单一 (depsrc, deprcv) 切片，维度变为 [1,1,north,east]""" + out = { + "depsrc": np.array([lib["depsrc"][isrc]], dtype=float), + "deprcv": np.array([lib["deprcv"][ircv]], dtype=float), + "north": lib["north"].copy(), + "east": lib["east"].copy(), + "src_va": np.array([lib["src_va"][isrc]], dtype=float), + "src_vb": np.array([lib["src_vb"][isrc]], dtype=float), + "src_rho": np.array([lib["src_rho"][isrc]], dtype=float), + "rcv_va": np.array([lib["rcv_va"][ircv]], dtype=float), + "rcv_vb": np.array([lib["rcv_vb"][ircv]], dtype=float), + "rcv_rho": np.array([lib["rcv_rho"][ircv]], dtype=float), + } + skip = set(out.keys()) | {"model"} + for k, v in lib.items(): + if k in skip: + continue + if getattr(v, "ndim", 0) != 4: + continue + # 通道数据 shape: (ndepsrc, ndeprcv, nnorth, neast) + out[k] = v[isrc:isrc + 1, ircv:ircv + 1, :, :].copy() + return out + + +def py_compute_to_nc(depsrcs, deprcvs, outpath: str): + pymod = pygrt.PyModel1D(MODNAME) + pymod.set_static_grn_path(outpath) + pymod.compute_static_grn( + depsrc=depsrcs, deprcv=deprcvs, norths=NORTHS, easts=EASTS, calc_upar=True, + ) + assert Path(outpath).is_file() + + +def load_syn_fields(path: Path) -> dict: + """读取合成 nc 中的场量（跳过坐标与介质）""" + out = {} + with netcdf_file(str(path), mmap=False) as f: + for k, v in f.variables.items(): + if k in SYN_SKIP: + continue + out[k] = np.array(v[:], dtype=float).copy() + return out + + +def compare_syn_fields(a: dict, b: dict, label: str) -> float: + """比较两个合成场量字典，返回平均相对误差""" + keys = sorted(set(a.keys()) & set(b.keys())) + if not keys: + raise ValueError(f"{label}: no syn fields to compare") + errors = [] + for k in keys: + va, vb = a[k], b[k] + if va.shape != vb.shape: + raise ValueError(f"{label}: shape mismatch on '{k}': {va.shape} vs {vb.shape}") + if np.all(va == 0.0) and np.all(vb == 0.0): + continue + denom = np.mean(np.abs(vb)) + if denom == 0.0: + denom = np.mean(np.abs(va)) + rerr = np.sum(np.abs(va - vb)) / denom + print(f" {k}: {rerr:.6e}") + errors.append(rerr) + mean_err = float(np.mean(errors)) if errors else 0.0 + print(f"{label}: mean relative error = {mean_err:.6e}") + if mean_err > MEAN_RERR_MAX: + raise ValueError(f"{label}: mean relative error too large ({mean_err})") + return mean_err + + +def c_static_syn(grn: Path, out: Path, extra: list) -> None: + run_grt(["static", "syn", f"-G{grn}", f"-S{format_float(SCALE)}", f"-O{out}", *extra, "-e"]) + + +def py_static_syn(grn: Path, out: Path, **kwargs) -> None: + model = pygrt.PyModel1D(MODNAME) + model.set_static_grn_path(grn) + model.compute_static_syn(scale=SCALE, output_path=out, calc_upar=True, **kwargs) + + +def compare_syn_cli_py() -> list: + """同一多深度库上，CLI 与 Python 合成结果应一致""" + errors = [] + grn = CMPDIR / "stgrn_mm.nc" + + print("\n--- syn exact depth Ds=2 Dr=0.5 ---") + c_out = CMPDIR / "stsyn_exact_c.nc" + py_out = CMPDIR / "stsyn_exact_py.nc" + c_static_syn(grn, c_out, ["-Ds2", "-Dr0.5"]) + py_static_syn(grn, py_out, depsrc=2.0, deprcv=0.5) + errors.append(compare_nc_files(py_out, c_out)) + + print("\n--- syn interpolated Ds=1.5 Dr=0.25 + new XY ---") + c_out = CMPDIR / "stsyn_interp_c.nc" + py_out = CMPDIR / "stsyn_interp_py.nc" + xy = ["-Ds1.5", "-Dr0.25", "-X-1/1/1", "-Y-1/1/1"] + c_static_syn(grn, c_out, xy) + py_static_syn( + grn, py_out, depsrc=1.5, deprcv=0.25, + norths=(-1.0, 1.0, 1.0), easts=(-1.0, 1.0, 1.0), + ) + errors.append(compare_nc_files(py_out, c_out)) + + print("\n--- syn -Q points ---") + rcv = CMPDIR / "rcv_pts.txt" + rcv.write_text("0 0 0\n1 1 0.5\n-1 0.5 1\n") + c_out = CMPDIR / "stsyn_q_c.nc" + py_out = CMPDIR / "stsyn_q_py.nc" + c_static_syn(grn, c_out, ["-Ds2", f"-Q{rcv}"]) + py_static_syn(grn, py_out, depsrc=2.0, recv_points=rcv) + errors.append(compare_nc_files(py_out, c_out)) + + print("\n--- syn + strain/rotation/stress ---") + c_ten = CMPDIR / "stsyn_ten_c.nc" + py_ten = CMPDIR / "stsyn_ten_py.nc" + c_static_syn(grn, c_ten, ["-Ds2", "-Dr0.5", "-N"]) + py_static_syn(grn, py_ten, depsrc=2.0, deprcv=0.5, zne=True) + run_grt(["static", "strain", str(c_ten)]) + run_grt(["static", "rotation", str(c_ten)]) + run_grt(["static", "stress", str(c_ten)]) + pygrt.utils.compute_strain(py_ten) + pygrt.utils.compute_rotation(py_ten) + pygrt.utils.compute_stress(py_ten) + errors.append(compare_nc_files(py_ten, c_ten)) + return errors + + +def compare_syn_multi_vs_single() -> list: + """多深度库在节点深度上的合成，应与对应单深度库一致""" + print("\n--- syn multi-depth lib at Ds=2 Dr=0.5 vs single-depth lib ---") + mm_out = CMPDIR / "stsyn_mm_node.nc" + ref_out = CMPDIR / "stsyn_ref_node.nc" + c_static_syn(CMPDIR / "stgrn_mm.nc", mm_out, ["-Ds2", "-Dr0.5"]) + c_static_syn(CMPDIR / "stgrn_ref_zs2_zr0p5.nc", ref_out, []) + return [compare_syn_fields(load_syn_fields(mm_out), load_syn_fields(ref_out), + "syn multi vs single [zs=2,zr=0.5]")] + + +def compare_syn_depth_linearity() -> list: + """线性插值：syn(1.5) 应等于 0.5*(syn(1)+syn(2))（同一 deprcv，不用 -Su）""" + print("\n--- syn depth interpolation linearity Ds=1.5 vs 0.5*(1+2) ---") + grn = CMPDIR / "stgrn_mm.nc" + out1 = CMPDIR / "stsyn_lin_z1.nc" + out2 = CMPDIR / "stsyn_lin_z2.nc" + outm = CMPDIR / "stsyn_lin_z15.nc" + c_static_syn(grn, out1, ["-Ds1", "-Dr0"]) + c_static_syn(grn, out2, ["-Ds2", "-Dr0"]) + c_static_syn(grn, outm, ["-Ds1.5", "-Dr0"]) + a = load_syn_fields(out1) + b = load_syn_fields(out2) + mid = load_syn_fields(outm) + pred = {k: 0.5 * (a[k] + b[k]) for k in mid if k in a and k in b} + return [compare_syn_fields(pred, mid, "syn depth linearity [Ds=1.5]")] + + +def main(): + all_errs = [] + + cases = [ + ("ss", 2.0, 0.0), + ("ms", DEPSRCS, 0.0), + ("mr", 2.0, DEPRCVS), + ("mm", DEPSRCS, DEPRCVS), + ] + + print("================ CLI vs Python ================") + for tag, zs, zr in cases: + c_path = CMPDIR / f"stgrn_{tag}.nc" + py_path = CMPDIR / f"stgrn_{tag}_py.nc" + print(f"\n--- case {tag} ---") + py_compute_to_nc(zs, zr, str(py_path)) + c_data = load_stgrnlib_nc(str(c_path)) + py_data = load_stgrnlib_nc(str(py_path)) + all_errs.append(compare_nc_data(c_data, py_data, f"CLI vs Python [{tag}]")) + + print("\n================ Multi vs Single (CLI) ================") + mm = load_stgrnlib_nc(str(CMPDIR / "stgrn_mm.nc")) + for isrc, zs in enumerate(DEPSRCS): + for ircv, zr in enumerate(DEPRCVS): + ref_path = CMPDIR / f"stgrn_ref_zs{_coord_tag(zs)}_zr{_coord_tag(zr)}.nc" + print(f"\n--- slice depsrc={zs}, deprcv={zr} ---") + sl = extract_slice(mm, isrc, ircv) + ref = load_stgrnlib_nc(str(ref_path)) + all_errs.append(compare_nc_data(sl, ref, f"multi vs single CLI [zs={zs},zr={zr}]")) + + print("\n================ Multi vs Single (Python) ================") + mm_py_path = CMPDIR / "stgrn_mm_py.nc" + mm_py = load_stgrnlib_nc(str(mm_py_path)) + for isrc, zs in enumerate(DEPSRCS): + for ircv, zr in enumerate(DEPRCVS): + single_path = CMPDIR / f"stgrn_py_zs{_coord_tag(zs)}_zr{_coord_tag(zr)}.nc" + print(f"\n--- py slice depsrc={zs}, deprcv={zr} ---") + py_compute_to_nc(float(zs), float(zr), str(single_path)) + sl = extract_slice(mm_py, isrc, ircv) + ref = load_stgrnlib_nc(str(single_path)) + all_errs.append(compare_nc_data(sl, ref, f"multi vs single Python [zs={zs},zr={zr}]")) + + print("\n================ Syn: CLI vs Python ================") + all_errs.extend(compare_syn_cli_py()) + + print("\n================ Syn: multi vs single (exact depth) ================") + all_errs.extend(compare_syn_multi_vs_single()) + + print("\n================ Syn: depth interpolation linearity ================") + all_errs.extend(compare_syn_depth_linearity()) + + all_errs = np.array(all_errs, dtype=float) + print("\n================ Summary ================") + print(all_errs) + print(f"mean={np.mean(all_errs):.6e} min={np.min(all_errs):.6e} max={np.max(all_errs):.6e}") + if np.mean(all_errs) > MEAN_RERR_MAX: + raise ValueError("overall mean relative error too large") + print("compare_stgrnlib.py: all checks passed") + + +if __name__ == "__main__": + main() diff --git a/test/_compare_c_py/test_cli_args.py b/test/_compare_c_py/test_cli_args.py index 2eb90350..34cf6932 100644 --- a/test/_compare_c_py/test_cli_args.py +++ b/test/_compare_c_py/test_cli_args.py @@ -7,6 +7,7 @@ from __future__ import annotations import sys +import warnings from pathlib import Path import pygrt @@ -284,6 +285,21 @@ def test_compute_static_grn_xy_and_distarr(): ) assert "-X" not in " ".join(cmd) assert "-Y" not in " ".join(cmd) + + # 多深度：应拼出 -Ds/-Dr，且 stats 被忽略 + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + model.compute_static_grn( + depsrc=[1.0, 2.0, 3.0], + deprcv=[0.0, 0.5], + norths=[-2.0, 2.0, 1.0], + easts=[-2.0, 2.0, 1.0], + stats=True, + ) + assert any("stats" in str(w.message) for w in caught) + cmd = runner.commands[-1] + assert_command_has(cmd, "static", "greenfn", "-Ds1,2,3", "-Dr0,0.5", "-X-2/2/1", "-Y-2/2/1") + assert not any(str(tok) == "-S" for tok in cmd) finally: _restore_run_grt(pygrt.pymod, original) @@ -421,8 +437,8 @@ def test_compute_static_syn_and_tensor_postprocess_args(): "static", "syn", f"-G{model.static_grn_path}", - "-Su1e+24", f"-O{out}", + "-Su1e+24", "-M33/50/120", "-X-5/5/1", "-Y-4/4/2", @@ -431,6 +447,59 @@ def test_compute_static_syn_and_tensor_postprocess_args(): ], ) + # 多深度点源 + 任意接收点 + model.compute_static_syn( + scale=1e20, + output_path=out, + source="EX", + depsrc=2.0, + recv_points=HERE / "rcv.txt", + ) + assert_command_has( + runner.commands[-1], + "static", "syn", + f"-G{model.static_grn_path}", + f"-O{out}", + "-S1e+20", + "-Ds2", + f"-Q{HERE / 'rcv.txt'}", + ) + + # 多深度点源 + 新网格 + 台站深度 + model.compute_static_syn( + scale=1e20, + output_path=out, + source="EX", + depsrc=2.0, + deprcv=0.5, + norths=[-2.0, 2.0, 1.0], + easts=[-2.0, 2.0, 1.0], + ) + assert_command_has( + runner.commands[-1], + "static", "syn", + "-Ds2", + "-Dr0.5", + "-X-2/2/1", + "-Y-2/2/1", + ) + + # 有限断层 + model.compute_static_syn( + output_path=out, + finite_fault=HERE / "cfaults.inp", + subfault_size=(1.0, 2.0), + calc_upar=True, + ) + assert_command_has( + runner.commands[-1], + "static", "syn", + f"-G{model.static_grn_path}", + f"-O{out}", + f"-C{HERE / 'cfaults.inp'}+i1/2", + "-e", + ) + # 张量后处理：目录走动态模块，文件走 static 模块 dyn = HERE / "_tmp_tensor_dyn" dyn.mkdir(exist_ok=True) diff --git a/test/_compare_c_py/test_compare_stgrnlib_c_py.sh b/test/_compare_c_py/test_compare_stgrnlib_c_py.sh new file mode 100644 index 00000000..375393cb --- /dev/null +++ b/test/_compare_c_py/test_compare_stgrnlib_c_py.sh @@ -0,0 +1,52 @@ +#!/bin/bash +# Compare multi/single-depth STGRNLIB: CLI vs Python, and multi vs single slices + +set -euo pipefail + +modname="milrow" + +# 粗网格以加快多深度测试 +x1=-2; x2=2; dx=1 +y1=-2; y2=2; dy=1 + +rm -rf stgrnlib_cmp +mkdir -p stgrnlib_cmp +cd stgrnlib_cmp + +# (1) 单震源 + 单台站 +grt static greenfn -M../../${modname} -D2/0 \ + -X$x1/$x2/$dx -Y$y1/$y2/$dy -e -Ostgrn_ss.nc + +# (2) 多震源 + 单台站 +grt static greenfn -M../../${modname} -Ds1,2,3 -Dr0 \ + -X$x1/$x2/$dx -Y$y1/$y2/$dy -e -Ostgrn_ms.nc + +# (3) 单震源 + 多台站 +grt static greenfn -M../../${modname} -Ds2 -Dr0,0.5,1 \ + -X$x1/$x2/$dx -Y$y1/$y2/$dy -e -Ostgrn_mr.nc + +# (4) 多震源 + 多台站 +grt static greenfn -M../../${modname} -Ds1,2,3 -Dr0,0.5,1 \ + -X$x1/$x2/$dx -Y$y1/$y2/$dy -e -Ostgrn_mm.nc + +# 深度标签：整数写 0/1，小数点换成 p（与 compare_stgrnlib.py 一致） +depth_tag() { + python -c "z=float('$1'); print(str(int(z)) if z==int(z) else str(z).replace('.','p'))" +} + +# 用于多 vs 单逐层对比的参考单深度结果 +for zs in 1 2 3; do + for zr in 0 0.5 1; do + zs_tag=$(depth_tag "$zs") + zr_tag=$(depth_tag "$zr") + grt static greenfn -M../../${modname} -D${zs}/${zr} \ + -X$x1/$x2/$dx -Y$y1/$y2/$dy -e \ + -Ostgrn_ref_zs${zs_tag}_zr${zr_tag}.nc + done +done + +cd - + +python -u compare_stgrnlib.py + +rm -rf stgrnlib_cmp diff --git a/test/static_greenfn/test_static_greenfn.py b/test/static_greenfn/test_static_greenfn.py index 2acf8cad..21e5cee4 100644 --- a/test/static_greenfn/test_static_greenfn.py +++ b/test/static_greenfn/test_static_greenfn.py @@ -1,7 +1,9 @@ import shutil +import warnings from pathlib import Path import pygrt +from scipy.io import netcdf_file depsrc = 2.0 deprcv = 0.0 @@ -41,6 +43,12 @@ depsrc=depsrc, deprcv=deprcv, norths=norths, easts=easts, k0=4, stats=True, ) +# 单深度输出应为 4D STGRNLIB（各深度维长度为 1） +with netcdf_file("stgrn.nc", mmap=False) as f: + assert "depsrc" in f.dimensions and "deprcv" in f.dimensions + assert f.dimensions["depsrc"] == 1 + assert f.dimensions["deprcv"] == 1 + # boundary condition pymod = pygrt.PyModel1D(modname, topbound="free", botbound="free") pymod.set_static_grn_path("stgrn.nc") @@ -54,9 +62,84 @@ pymod.set_static_grn_path("stgrn.nc") pymod.compute_static_grn(depsrc=depsrc, deprcv=deprcv, norths=norths, easts=easts) -for name in ["stgrn.nc", "stgrtstats"]: +# -------------------- 多深度功能 -------------------- +norths_c = [-2.0, 2.0, 1.0] +easts_c = [-2.0, 2.0, 1.0] +depsrcs = [1.0, 2.0, 3.0] +deprcvs = [0.0, 0.5] + +pymod_m = pygrt.PyModel1D(modname) +pymod_m.set_static_grn_path("stgrn_py_multi.nc") +pymod_m.compute_static_grn(depsrc=depsrcs, deprcv=deprcvs, norths=norths_c, easts=easts_c) +assert Path("stgrn_py_multi.nc").is_file() +with netcdf_file("stgrn_py_multi.nc", mmap=False) as f: + assert f.dimensions["depsrc"] == 3 + assert f.dimensions["deprcv"] == 2 + +# 仅多震源深度 +pymod_ms = pygrt.PyModel1D(modname) +pymod_ms.set_static_grn_path("stgrn_py_ms.nc") +pymod_ms.compute_static_grn( + depsrc=[1.0, 2.0], deprcv=0.0, norths=norths_c, easts=easts_c, calc_upar=True, +) + +# 仅多台站深度 +pymod_mr = pygrt.PyModel1D(modname) +pymod_mr.set_static_grn_path("stgrn_py_mr.nc") +pymod_mr.compute_static_grn(depsrc=2.0, deprcv=[0.0, 0.5], norths=norths_c, easts=easts_c) + +# -R / distarr 建库 +pymod_r = pygrt.PyModel1D(modname) +pymod_r.set_static_grn_path("stgrn_py_r.nc") +pymod_r.compute_static_grn( + depsrc=2.0, deprcv=0.0, distarr=[0.0, 1.0, 2.0, 4.0], +) +with netcdf_file("stgrn_py_r.nc", mmap=False) as f: + assert f.dimensions["north"] == 1 + assert f.dimensions["east"] == 4 + +# -------------------- 错误 / 警告（Python）-------------------- +try: + pymod_m.compute_static_grn(depsrc=-1.0, deprcv=0.0, norths=norths_c, easts=easts_c) + raise AssertionError("negative depsrc should raise") +except ValueError: + pass + +try: + pymod_m.compute_static_grn(depsrc=1.0, deprcv=-0.5, norths=norths_c, easts=easts_c) + raise AssertionError("negative deprcv should raise") +except ValueError: + pass + +try: + pymod_m.compute_static_grn(depsrc=[], deprcv=0.0, norths=norths_c, easts=easts_c) + raise AssertionError("empty depsrc should raise") +except ValueError: + pass + +try: + pymod_r.compute_static_grn(depsrc=2.0, deprcv=0.0, distarr=[0.0, 2.0, 1.0]) + raise AssertionError("non-ascending distarr should raise") +except ValueError: + pass + +# 多深度 stats 应警告并忽略 +with warnings.catch_warnings(record=True) as w: + warnings.simplefilter("always") + pymod_m.set_static_grn_path("stgrn_py_multi2.nc") + pymod_m.compute_static_grn( + depsrc=depsrcs, deprcv=deprcvs, norths=norths_c, easts=easts_c, stats=True, + ) + assert any("stats" in str(x.message) for x in w) + +for name in [ + "stgrn.nc", "stgrn_py_multi.nc", "stgrn_py_ms.nc", "stgrn_py_mr.nc", + "stgrn_py_multi2.nc", "stgrn_py_r.nc", "stgrtstats", +]: p = Path(name) if p.is_dir(): shutil.rmtree(p, ignore_errors=True) elif p.is_file(): p.unlink(missing_ok=True) + +print("test_static_greenfn.py: all checks passed") diff --git a/test/static_greenfn/test_static_greenfn.sh b/test/static_greenfn/test_static_greenfn.sh index 6e231088..4c59eb2c 100755 --- a/test/static_greenfn/test_static_greenfn.sh +++ b/test/static_greenfn/test_static_greenfn.sh @@ -2,9 +2,52 @@ set -euo pipefail +expect_fail() { + # 期望命令失败；成功则报错退出 + local desc="$1" + shift + set +e + "$@" >/dev/null 2>&1 + local ret=$? + set -e + if [ "$ret" -eq 0 ]; then + echo "ERROR: expected failure but succeeded: $desc" >&2 + exit 1 + fi + echo "OK (failed as expected): $desc" +} + +expect_warn() { + # 期望命令成功且 stderr 含指定关键字 + local desc="$1" + local key="$2" + shift 2 + local tmp + tmp=$(mktemp) + set +e + "$@" >"$tmp" 2>&1 + local ret=$? + set -e + if [ "$ret" -ne 0 ]; then + echo "ERROR: expected success with warning but failed: $desc" >&2 + cat "$tmp" >&2 + rm -f "$tmp" + exit 1 + fi + if ! grep -q "$key" "$tmp"; then + echo "ERROR: expected warning containing '$key': $desc" >&2 + cat "$tmp" >&2 + rm -f "$tmp" + exit 1 + fi + rm -f "$tmp" + echo "OK (warned as expected): $desc" +} + grt static greenfn -h grt static_greenfn -h +# -------------------- 单深度（旧 -D）-------------------- grt static greenfn -M../milrow -D2/0 -X-3/3/0.2 -Y-2/2/0.2 -Ostgrn.nc grt static greenfn -M../milrow -D2/0 -X-3/3/0.2 -Y-2/2/0.2 -e -Ostgrn.nc grt static greenfn -M../milrow -D2/0 -X-3/3/0.2 -Y-2/2/0.2 -L20 -Ostgrn.nc @@ -23,12 +66,59 @@ grt static greenfn -M../milrow -D2/0 -Rdists -Ostgrn.nc rm -rf dists grt static greenfn -M../milrow -D2/0 -R2,3,5,8 -Ostgrn.nc - # boundary grt static greenfn -M../milrow -D2/0 -X-3/3/0.2 -Y-2/2/0.2 -BrF -Ostgrn.nc grt static greenfn -M../milrow -D2/0 -X-3/3/0.2 -Y-2/2/0.2 -BhR -Ostgrn.nc grt static greenfn -M../milrow -D2/0 -X-3/3/0.2 -Y-2/2/0.2 -BrH -Ostgrn.nc +# -------------------- 多深度（-Ds / -Dr）-------------------- +# 逗号列表 +grt static greenfn -M../milrow -Ds1,2,3 -Dr0 -X-2/2/1 -Y-2/2/1 -Ostgrn_multi.nc +# 等间距 +grt static greenfn -M../milrow -Ds1/3/1 -Dr0/1/1 -X-2/2/1 -Y-2/2/1 -Ostgrn_multi.nc +# 单深度也可用 -Ds/-Dr +grt static greenfn -M../milrow -Ds2 -Dr0 -X-2/2/1 -Y-2/2/1 -Ostgrn_multi.nc +# 带位移偏导 +grt static greenfn -M../milrow -Ds1,2 -Dr0,0.5 -X-2/2/1 -Y-2/2/1 -e -Ostgrn_multi.nc +# 从文件读深度 +printf "1\n2\n3\n" > depsrc_list +printf "0\n0.5\n" > deprcv_list +grt static greenfn -M../milrow -Dsdepsrc_list -Drdeprcv_list -X-2/2/1 -Y-2/2/1 -Ostgrn_multi.nc +rm -f depsrc_list deprcv_list + +# 多深度时 -S 应警告并忽略，但仍成功 +expect_warn "multi-depth -S ignored" "ignored" \ + grt static greenfn -M../milrow -Ds1,2 -Dr0 -X-2/2/1 -Y-2/2/1 -S -Ostgrn_multi.nc + +# -------------------- 错误参数（CLI）-------------------- +expect_fail "-D and -Ds/-Dr mutually exclusive" \ + grt static greenfn -M../milrow -D2/0 -Ds1,2 -Dr0 -X-2/2/1 -Y-2/2/1 -Ostgrn_bad.nc + +expect_fail "-Ds without -Dr" \ + grt static greenfn -M../milrow -Ds1,2 -X-2/2/1 -Y-2/2/1 -Ostgrn_bad.nc + +expect_fail "-Dr without -Ds" \ + grt static greenfn -M../milrow -Dr0 -X-2/2/1 -Y-2/2/1 -Ostgrn_bad.nc + +expect_fail "missing depth option" \ + grt static greenfn -M../milrow -X-2/2/1 -Y-2/2/1 -Ostgrn_bad.nc + +expect_fail "negative depth in -Ds" \ + grt static greenfn -M../milrow -Ds-1,2 -Dr0 -X-2/2/1 -Y-2/2/1 -Ostgrn_bad.nc + +expect_fail "negative depth in -D" \ + grt static greenfn -M../milrow -D-1/0 -X-2/2/1 -Y-2/2/1 -Ostgrn_bad.nc + +expect_fail "nonpositive depth spacing" \ + grt static greenfn -M../milrow -Ds1/3/0 -Dr0 -X-2/2/1 -Y-2/2/1 -Ostgrn_bad.nc + +expect_fail "depth start > end" \ + grt static greenfn -M../milrow -Ds3/1/1 -Dr0 -X-2/2/1 -Y-2/2/1 -Ostgrn_bad.nc + +expect_fail "non-ascending -R list" \ + grt static greenfn -M../milrow -D2/0 -R3,1,2 -Ostgrn_bad.nc + +# -------------------- Python -------------------- python -u test_static_greenfn.py -rm -rf stgrt* *.nc \ No newline at end of file +rm -rf stgrt* *.nc diff --git a/test/static_syn/test_static_syn.py b/test/static_syn/test_static_syn.py index 91393c27..2534c695 100644 --- a/test/static_syn/test_static_syn.py +++ b/test/static_syn/test_static_syn.py @@ -1,6 +1,8 @@ +import shutil from pathlib import Path import pygrt +from scipy.io import netcdf_file depsrc = 2.0 deprcv = 0.0 @@ -32,7 +34,136 @@ moment_tensor=(1, -2, -5, 0.5, 3, 1.2), ) -for name in ["stgrn.nc", "stsyn.nc"]: +# ZNE / 空间导数 / 新 XY 网格 +pymod.compute_static_syn( + scale=1e20, output_path="stsyn.nc", source="EX", zne=True, calc_upar=True, +) +pymod.compute_static_syn( + scale=1e20, output_path="stsyn_xy.nc", source="EX", + norths=[-2.0, 2.0, 0.5], easts=[-1.0, 1.0, 0.5], +) + +# -------------------- -R 建库后合成 -------------------- +pymod_r = pygrt.PyModel1D(modname) +pymod_r.set_static_grn_path("stgrn_r.nc") +pymod_r.compute_static_grn( + depsrc=depsrc, deprcv=deprcv, distarr=[0.0, 1.0, 2.0, 4.0, 8.0], calc_upar=True, +) +pymod_r.compute_static_syn(scale=1e20, output_path="stsyn_r.nc", source="EX") +pymod_r.compute_static_syn( + scale=1e20, output_path="stsyn_rxy.nc", source="EX", + norths=[-3.0, 3.0, 1.0], easts=[-2.0, 2.0, 1.0], calc_upar=True, +) + +# -------------------- 多深度：点源 -Ds、深度插值、-Q、有限断层 -------------------- +pymod_m = pygrt.PyModel1D(modname) +pymod_m.set_static_grn_path("stgrn_md.nc") +pymod_m.compute_static_grn( + depsrc=[1.0, 2.0, 3.0], deprcv=0.0, + distarr=[0.0, 1.0, 2.0, 4.0, 8.0], calc_upar=True, +) +pymod_m.compute_static_syn( + scale=1e16, output_path="stsyn_md.nc", source="EX", depsrc=2.0, scale_with_mu=True, +) +pymod_m.compute_static_syn( + scale=1e16, output_path="stsyn_interp.nc", source="EX", depsrc=1.5, + norths=[-2.0, 2.0, 1.0], easts=[-2.0, 2.0, 1.0], + scale_with_mu=True, calc_upar=True, +) + +rcv = Path("rcv_pts.txt") +rcv.write_text("# north east depth (km)\n0 0 0\n1 2 0\n-1 1 0\n") +pymod_m.compute_static_syn( + scale=1e16, output_path="stsyn_q.nc", source="EX", depsrc=2.0, recv_points=rcv, +) +with netcdf_file("stsyn_q.nc", mmap=False) as f: + assert "point" in f.dimensions + assert f.dimensions["point"] == 3 + +ff = Path("cfaults_tiny.inp") +# W=(2.8-1.2)/sin(90°)=1.6 km，dW=1 → 末块短于 dW，覆盖余数子断层中心 +ff.write_text( + " # X-start Y-start X-fin Y-fin Kode shear(m) reverse(m) dip angle top(km) bot(km)\n" + "xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx\n" + " 1 0.0000 0.0000 2.0000 0.0000 100 0.1000 0.0000 90.00 1.2000 2.8000\n" +) +pymod_m.compute_static_syn( + output_path="stsyn_ff.nc", finite_fault=ff, subfault_size=(1.0, 1.0), +) + +# 多台站深度必须给 deprcv +pymod_mr = pygrt.PyModel1D(modname) +pymod_mr.set_static_grn_path("stgrn_mr.nc") +pymod_mr.compute_static_grn( + depsrc=2.0, deprcv=[0.0, 0.5], distarr=[0.0, 5.0], calc_upar=True, +) +pymod_mr.compute_static_syn( + scale=1e20, output_path="stsyn_dr.nc", source="EX", + deprcv=0.25, norths=[-2.0, 2.0, 1.0], easts=[-2.0, 2.0, 1.0], +) + +# -------------------- 错误参数 -------------------- +try: + pymod_m.compute_static_syn( + scale=1e20, output_path="stsyn_bad.nc", source="EX", + recv_points=rcv, norths=[-1.0, 1.0, 1.0], easts=[-1.0, 1.0, 1.0], + ) + raise AssertionError("recv_points with norths/easts should raise") +except ValueError: + pass + +try: + pymod_m.compute_static_syn( + scale=1e20, output_path="stsyn_bad.nc", source="EX", + recv_points=rcv, deprcv=0.0, + ) + raise AssertionError("recv_points with deprcv should raise") +except ValueError: + pass + +try: + pymod_m.compute_static_syn( + scale=1e20, output_path="stsyn_bad.nc", source="EX", + finite_fault=ff, + ) + raise AssertionError("finite_fault with scale should raise") +except ValueError: + pass + +bad_dip = Path("cfaults_bad_dip.inp") +bad_dip.write_text( + " # X-start Y-start X-fin Y-fin Kode shear(m) reverse(m) dip angle top(km) bot(km)\n" + "xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx\n" + " 1 0.0000 0.0000 2.0000 0.0000 100 0.1000 0.0000 0.00 1.2000 2.8000\n" +) +try: + pymod_m.compute_static_syn(output_path="stsyn_bad.nc", finite_fault=bad_dip) + raise AssertionError("dip=0 should fail in C") +except RuntimeError: + pass + +bad_bot = Path("cfaults_bad_bot.inp") +bad_bot.write_text( + " # X-start Y-start X-fin Y-fin Kode shear(m) reverse(m) dip angle top(km) bot(km)\n" + "xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx\n" + " 1 0.0000 0.0000 2.0000 0.0000 100 0.1000 0.0000 90.00 2.8000 1.2000\n" +) +try: + pymod_m.compute_static_syn(output_path="stsyn_bad.nc", finite_fault=bad_bot) + raise AssertionError("bot < top should fail in C") +except RuntimeError: + pass + +for name in [ + "stgrn.nc", "stgrn_r.nc", "stgrn_md.nc", "stgrn_mr.nc", + "stsyn.nc", "stsyn_xy.nc", "stsyn_r.nc", "stsyn_rxy.nc", + "stsyn_md.nc", "stsyn_interp.nc", "stsyn_q.nc", "stsyn_ff.nc", "stsyn_dr.nc", + "rcv_pts.txt", "cfaults_tiny.inp", "cfaults_bad_dip.inp", "cfaults_bad_bot.inp", +]: p = Path(name) - if p.is_file(): + if p.is_dir(): + shutil.rmtree(p, ignore_errors=True) + elif p.is_file(): p.unlink(missing_ok=True) + +print("test_static_syn.py: all checks passed") diff --git a/test/static_syn/test_static_syn.sh b/test/static_syn/test_static_syn.sh index c9d37044..2b053ef8 100755 --- a/test/static_syn/test_static_syn.sh +++ b/test/static_syn/test_static_syn.sh @@ -2,9 +2,25 @@ set -euo pipefail +expect_fail() { + # 期望命令失败；成功则报错退出 + local desc="$1" + shift + set +e + "$@" >/dev/null 2>&1 + local ret=$? + set -e + if [ "$ret" -eq 0 ]; then + echo "ERROR: expected failure but succeeded: $desc" >&2 + exit 1 + fi + echo "OK (failed as expected): $desc" +} + grt static syn -h grt static_syn -h +# -------------------- 单深度 XY 库 -------------------- grt static greenfn -M../milrow -D2/0 -X-3/3/0.2 -Y-2/2/0.2 -e -Ostgrn.nc grt static syn -S1e20 -Gstgrn.nc -Ostsyn.nc @@ -17,6 +33,82 @@ grt static syn -S1e20 -T1/-2/-5/0.5/3/1.2 -Gstgrn.nc -Ostsyn.nc grt static syn -S1e20 -F2/-1/4 -e -Gstgrn.nc -Ostsyn.nc grt static syn -S1e20 -F2/-1/4 -N -e -Gstgrn.nc -Ostsyn.nc +# 合成时指定新 XY 网格（震中距插值） +grt static syn -S1e20 -Gstgrn.nc -X-2/2/0.5 -Y-1/1/0.5 -Ostsyn_xy.nc + +# -------------------- -R 建库后合成 -------------------- +grt static greenfn -M../milrow -D2/0 -R0/8/1 -e -Ostgrn_r.nc +# 延用库的震中距轴 +grt static syn -S1e20 -Gstgrn_r.nc -Ostsyn_r.nc +# 指定新 XY 网格 +grt static syn -S1e20 -Gstgrn_r.nc -X-3/3/1 -Y-2/2/1 -e -Ostsyn_rxy.nc + +# -------------------- 多深度库：点源 -Ds；深度插值；-Q -------------------- +grt static greenfn -M../milrow -Ds1,2,3 -Dr0 -R0/8/1 -e -Ostgrn_md.nc +grt static syn -Gstgrn_md.nc -Su1e16 -Ds2 -Ostsyn_md.nc +# 震源深度插值（库节点之间） +grt static syn -Gstgrn_md.nc -Su1e16 -Ds1.5 -X-2/2/1 -Y-2/2/1 -e -Ostsyn_interp.nc + +cat > rcv_pts.txt <<'EOF' +# north east depth (km) +0 0 0 +1 2 0 +-1 1 0 +EOF +grt static syn -Gstgrn_md.nc -Su1e16 -Ds2 -Qrcv_pts.txt -Ostsyn_q.nc + +# 多台站深度：必须 -Dr +grt static greenfn -M../milrow -Ds2 -Dr0,0.5 -R0,5 -e -Ostgrn_mr.nc +grt static syn -Gstgrn_mr.nc -S1e20 -Dr0.25 -X-2/2/1 -Y-2/2/1 -Ostsyn_dr.nc + +# -------------------- 有限断层（库震源深度覆盖断层 top/bot）-------------------- +# W=(2.8-1.2)/sin(90°)=1.6 km，dW=1 → 末块短于 dW，用于覆盖余数子断层中心 +cat > cfaults_tiny.inp <<'EOF' + # X-start Y-start X-fin Y-fin Kode shear(m) reverse(m) dip angle top(km) bot(km) +xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx + 1 0.0000 0.0000 2.0000 0.0000 100 0.1000 0.0000 90.00 1.2000 2.8000 +EOF +grt static syn -Gstgrn_md.nc -Ccfaults_tiny.inp+i1/1 -Ostsyn_ff.nc +grt static syn -Gstgrn_md.nc -Ccfaults_tiny.inp+i1/1 -e -Qrcv_pts.txt -Ostsyn_ffq.nc + +# -------------------- 错误参数 -------------------- +expect_fail "multi-src library requires -Ds" \ + grt static syn -Gstgrn_md.nc -S1e20 -Ostsyn_bad.nc + +expect_fail "multi-rcv library requires -Dr" \ + grt static syn -Gstgrn_mr.nc -S1e20 -Ostsyn_bad.nc + +expect_fail "single-rcv library forbids -Dr" \ + grt static syn -Gstgrn.nc -S1e20 -Dr0 -Ostsyn_bad.nc + +expect_fail "-Q mutually exclusive with -X/-Y" \ + grt static syn -Gstgrn_md.nc -S1e20 -Ds2 -Qrcv_pts.txt -X-1/1/1 -Y-1/1/1 -Ostsyn_bad.nc + +expect_fail "-Q mutually exclusive with -Dr" \ + grt static syn -Gstgrn_md.nc -S1e20 -Ds2 -Dr0 -Qrcv_pts.txt -Ostsyn_bad.nc + +expect_fail "finite fault requires ndepsrc>1" \ + grt static syn -Gstgrn.nc -Ccfaults_tiny.inp -Ostsyn_bad.nc + +expect_fail "finite fault mutually exclusive with -S" \ + grt static syn -Gstgrn_md.nc -S1e20 -Ccfaults_tiny.inp -Ostsyn_bad.nc + +cat > cfaults_bad_dip.inp <<'EOF' + # X-start Y-start X-fin Y-fin Kode shear(m) reverse(m) dip angle top(km) bot(km) +xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx + 1 0.0000 0.0000 2.0000 0.0000 100 0.1000 0.0000 0.00 1.2000 2.8000 +EOF +expect_fail "finite fault dip must be in (0, 90]" \ + grt static syn -Gstgrn_md.nc -Ccfaults_bad_dip.inp -Ostsyn_bad.nc + +cat > cfaults_bad_bot.inp <<'EOF' + # X-start Y-start X-fin Y-fin Kode shear(m) reverse(m) dip angle top(km) bot(km) +xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx xxxxxxxxxx + 1 0.0000 0.0000 2.0000 0.0000 100 0.1000 0.0000 90.00 2.8000 1.2000 +EOF +expect_fail "finite fault bot must be greater than top" \ + grt static syn -Gstgrn_md.nc -Ccfaults_bad_bot.inp -Ostsyn_bad.nc + python -u test_static_syn.py -rm -rf *.nc \ No newline at end of file +rm -rf *.nc rcv_pts.txt cfaults_tiny.inp cfaults_bad_dip.inp cfaults_bad_bot.inp diff --git a/test/tensors/test_tensors.py b/test/tensors/test_tensors.py index 21058f45..ee11ceb8 100644 --- a/test/tensors/test_tensors.py +++ b/test/tensors/test_tensors.py @@ -33,7 +33,40 @@ pygrt.utils.compute_stress("syn_zne") pygrt.utils.compute_rotation("syn_zne") -for name in ["GRN", "syn", "syn_zne"]: +# -------------------- 静态应变 / 应力 / 旋转 -------------------- +pymod_s = pygrt.PyModel1D(modname) +pymod_s.set_static_grn_path("stgrn.nc") +pymod_s.compute_static_grn( + depsrc=2.0, deprcv=0.0, norths=[-3.0, 3.0, 1.0], easts=[-2.0, 2.0, 1.0], + calc_upar=True, +) +pymod_s.compute_static_syn( + scale=1e20, output_path="stsyn.nc", source="EX", calc_upar=True, +) +pygrt.utils.compute_strain("stsyn.nc") +pygrt.utils.compute_stress("stsyn.nc") +pygrt.utils.compute_rotation("stsyn.nc") + +pymod_s.compute_static_syn( + scale=1e20, output_path="stsyn_zne.nc", source="EX", zne=True, calc_upar=True, +) +pygrt.utils.compute_strain("stsyn_zne.nc") +pygrt.utils.compute_stress("stsyn_zne.nc") +pygrt.utils.compute_rotation("stsyn_zne.nc") + +rcv = Path("rcv_pts.txt") +rcv.write_text("# north east depth (km)\n0 0 0\n1 2 0\n-1 1 0\n") +pymod_s.compute_static_syn( + scale=1e20, output_path="stsyn_q.nc", source="EX", + recv_points=rcv, calc_upar=True, +) +pygrt.utils.compute_strain("stsyn_q.nc") +pygrt.utils.compute_stress("stsyn_q.nc") +pygrt.utils.compute_rotation("stsyn_q.nc") + +for name in ["GRN", "syn", "syn_zne", "stgrn.nc", "stsyn.nc", "stsyn_zne.nc", "stsyn_q.nc", "rcv_pts.txt"]: p = Path(name) if p.is_dir(): shutil.rmtree(p, ignore_errors=True) + elif p.is_file(): + p.unlink(missing_ok=True) diff --git a/test/tensors/test_tensors.sh b/test/tensors/test_tensors.sh index ec8cfa65..82258cbf 100755 --- a/test/tensors/test_tensors.sh +++ b/test/tensors/test_tensors.sh @@ -17,7 +17,34 @@ grt strain syn_ZNE grt stress syn_ZNE grt rotation syn_ZNE -python -u test_tensors.py +# -------------------- 静态应变 / 应力 / 旋转 -------------------- +grt static strain -h +grt static stress -h +grt static rotation -h + +grt static greenfn -M../milrow -D2/0 -X-3/3/1 -Y-2/2/1 -e -Ostgrn.nc +grt static syn -Gstgrn.nc -S1e20 -e -Ostsyn.nc +grt static strain stsyn.nc +grt static stress stsyn.nc +grt static rotation stsyn.nc + +grt static syn -Gstgrn.nc -S1e20 -e -N -Ostsyn_ZNE.nc +grt static strain stsyn_ZNE.nc +grt static stress stsyn_ZNE.nc +grt static rotation stsyn_ZNE.nc +# 任意接收点布局 +cat > rcv_pts.txt <<'EOF' +# north east depth (km) +0 0 0 +1 2 0 +-1 1 0 +EOF +grt static syn -Gstgrn.nc -S1e20 -e -Qrcv_pts.txt -Ostsyn_q.nc +grt static strain stsyn_q.nc +grt static stress stsyn_q.nc +grt static rotation stsyn_q.nc + +python -u test_tensors.py -rm -rf GRN syn* +rm -rf GRN syn* stgrn.nc stsyn*.nc rcv_pts.txt