diff --git a/docs/source/Tutorial/static/run/run.py b/docs/source/Tutorial/static/run/run.py index 09558d8..c6bf5df 100755 --- a/docs/source/Tutorial/static/run/run.py +++ b/docs/source/Tutorial/static/run/run.py @@ -9,7 +9,10 @@ xarr = np.linspace(-3, 3, 41) yarr = np.linspace(-2.5, 2.5, 33) -static_grn = pymod.compute_static_grn(xarr, yarr) +# 可以设置 distarr 来指定震中距序列 +# static_grn = pymod.compute_static_grn(distarr=np.arange(0,10+1e-8,0.1)) +# 也可以设置 xarr 和 yarr 来指定 XY 网格 +static_grn = pymod.compute_static_grn(xarr=xarr, yarr=yarr) print(static_grn.keys()) # dict_keys(['_xarr', '_yarr', '_src_va', '_src_vb', '_src_rho', '_rcv_va', '_rcv_vb', '_rcv_rho', 'EXZ', 'VFZ', 'DDZ', 'HFZ', 'DSZ', 'SSZ', 'EXR', 'VFR', 'DDR', 'HFR', 'DSR', 'SSR', 'HFT', 'DST', 'SST']) # END GRN @@ -114,3 +117,15 @@ def plot_static(static_syn:dict, out:Union[str,None]=None): plot_static(static_syn, "syn_mt2.svg") # END SYN MT2 # --------------------------------------------------------------------------------- + + +# --------------------------------------------------------------------------------- +# BEGIN NEW XY +xarr2 = np.arange(-3, 3+1e-8, 0.2) +yarr2 = np.arange(-2.5, 2.5+1e-8, 0.25) +static_syn = pygrt.utils.gen_syn_from_gf_DC(static_grn, M0=1e24, strike=33, dip=90, rake=0, ZNE=True, xarr=xarr2, yarr=yarr2) +print(static_syn.keys()) +# dict_keys(['_xarr', '_yarr', '_src_va', '_src_vb', '_src_rho', '_rcv_va', '_rcv_vb', '_rcv_rho', 'Z', 'N', 'E']) +plot_static(static_syn, "synXY_dc2.svg") +# END NEW XY +# --------------------------------------------------------------------------------- diff --git a/docs/source/Tutorial/static/run/run.sh b/docs/source/Tutorial/static/run/run.sh index c9bb59c..7b90264 100755 --- a/docs/source/Tutorial/static/run/run.sh +++ b/docs/source/Tutorial/static/run/run.sh @@ -41,6 +41,9 @@ y1=-2.5 y2=2.5 dy=0.15 # 输出到NetCDF网格文件 +# 可以使用 -R 来指定震中距序列 +# grt static greenfn -Mmilrow -D${depsrc}/${deprcv} -R0/10/0.1 -Ostgrn.nc +# 也可以使用 -X 和 -Y 来指定 XY 网格 grt static greenfn -Mmilrow -D${depsrc}/${deprcv} -X$x1/$x2/$dx -Y$y1/$y2/$dy -Ostgrn.nc # END GRN # --------------------------------------------------------------------------------- @@ -160,7 +163,6 @@ grt static syn -S1e24 -T0/-0.2/0/0/0/0 -N -Gstgrn.nc -Ostsyn_mt2.nc # END SYN MT2 # --------------------------------------------------------------------------------- - # X Y depth mrr mtt mff mrt mrf mtf exp # mzz mxx myy mxz -myz -mxy gmt begin syn_mt2 pdf @@ -171,6 +173,23 @@ EOF gmt colorbar -Bx+l"Z (cm)" gmt end + +# --------------------------------------------------------------------------------- +# BEGIN NEW XY +# 从网格文件中读取格林函数，再将合成结果写入新网格 +grt static syn -S1e24 -M33/90/0 -N -Gstgrn.nc -OstsynXY_dc2.nc -X-3/3/0.2 -Y-2.5/2.5/0.25 +# END NEW XY +# --------------------------------------------------------------------------------- + +gmt begin synXY_dc2 pdf + gmtplot_static stsynXY_dc2.nc -Si0.03c + gmt meca -Sa0.5c < -D/ -X// \n" -" -Y// -O [-L] [-C[d|p|n]] \n" -" [-Bf|F|r|R|h|H] [-K[+k][+e]] [-S] [-e]\n" +" grt static greenfn -M -D/ -O \n" +" [-X//] [-Y//] [-R//]\n" +" [-L] [-C[d|p|n]] [-Bf|F|r|R|h|H] \n" +" [-K[+k][+e]] [-S] [-e]\n" +"\n" +" There're two ways to define the \"epicentral distances\":\n" +" 1. set both -X and -Y to define a XY grid in advance.\n" +" 2. simply set -R, which equal to \"-X0/0/1 -Y//\".\n" "\n\n" "Options:\n" "----------------------------------------------------------------\n" @@ -184,6 +189,9 @@ printf("\n" " : end coordinate (km).\n" " : sampling interval (km).\n" "\n" +" -R//\n" +" equal to \"-X0/0/1 -Y//\". \n" +"\n" " -O Filepath to output nc grid.\n" "\n" " -L[][+l][+a][+o]\n" @@ -241,7 +249,8 @@ printf("\n" "\n\n" "Examples:\n" "----------------------------------------------------------------\n" -" grt static greenfn -Mmilrow -D2/0 -X-10/10/20 -Y-10/10/20 -Ostgrn.nc\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" "\n\n\n" ); } @@ -260,7 +269,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ int opt; - while ((opt = getopt(argc, argv, ":M:D:B:L:C:K:X:Y:O:Seh")) != -1) { + while ((opt = getopt(argc, argv, ":M:D:B:L:C:K:X:Y:R:O:Seh")) != -1) { switch (opt) { // 模型路径，其中每行分别为 // 厚度(km) Vp(km/s) Vs(km/s) Rho(g/cm^3) Qp Qs @@ -466,6 +475,37 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } break; + // -R 算是别名，相当于 -X0/0/1, -Yy1/y2/dy, 以此方式指定一维震中距序列 + case 'R': + Ctrl->X.active = Ctrl->Y.active = true; + { + real_t a1, a2, delta; + if(3 != sscanf(optarg, "%lf/%lf/%lf", &a1, &a2, &delta)){ + GRTBadOptionError(R, ""); + }; + if(delta <= 0){ + GRTBadOptionError(R, "Can't set nonpositive dr(%f)", delta); + } + // 特别检查 ys 不能是负数 + if(a1 < 0.0){ + GRTBadOptionError(R, "Can't set nonpositive r1(%f)", a1); + } + if(a1 > a2){ + GRTBadOptionError(R, "r1(%f) > r2(%f).", a1, a2); + } + + Ctrl->Y.ny = floor((a2-a1)/delta) + 1; + Ctrl->Y.ys = (real_t*)calloc(Ctrl->Y.ny, sizeof(real_t)); + for(size_t i=0; iY.ny; ++i){ + Ctrl->Y.ys[i] = a1 + delta*i; + } + + Ctrl->X.nx = 1; + Ctrl->X.xs = (real_t*)calloc(Ctrl->X.nx, sizeof(real_t)); + Ctrl->X.xs[0] = 0.0; + } + break; + // 输出 nc 文件名 case 'O': Ctrl->O.active = true; diff --git a/pygrt/C_extension/src/static/grt_static_syn.c b/pygrt/C_extension/src/static/grt_static_syn.c index 66fd351..91d7a82 100644 --- a/pygrt/C_extension/src/static/grt_static_syn.c +++ b/pygrt/C_extension/src/static/grt_static_syn.c @@ -12,6 +12,7 @@ #include "grt/common/radiation.h" #include "grt/common/coord.h" #include "grt/common/util.h" +#include "grt/common/search.h" #include "grt/common/mynetcdf.h" #include "grt.h" @@ -46,6 +47,18 @@ typedef struct { struct { bool active; } T; + /** X 坐标 */ + struct { + bool active; + size_t nx; + real_t *xs; + } X; + /** Y 坐标 */ + struct { + bool active; + size_t ny; + real_t *ys; + } Y; /** 输出 nc 文件名 */ struct { bool active; @@ -60,12 +73,12 @@ typedef struct { bool active; } e; + // 是否使用新网格 + bool isnewXYgrid; + // 存储不同震源的震源机制相关参数的数组 real_t mchn[GRT_MECHANISM_NUM]; - // 方向因子数组 - realChnlGrid srcRadi; - // 最终要计算的震源类型 GRT_SYN_TYPE computeType; char s_computeType[3]; @@ -77,6 +90,12 @@ static void free_Ctrl(GRT_MODULE_CTRL *Ctrl){ // G GRT_SAFE_FREE_PTR(Ctrl->G.s_ingrid); + // X + GRT_SAFE_FREE_PTR(Ctrl->X.xs); + + // Y + GRT_SAFE_FREE_PTR(Ctrl->Y.ys); + // O GRT_SAFE_FREE_PTR(Ctrl->O.s_outgrid); @@ -94,15 +113,18 @@ printf("\n" " + Radial Outward (R),\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 and -Y to define a new XY grid, the program \n" +" will automatically find the nearest Green's functions for each node.\n" "\n\n" "Usage:\n" "----------------------------------------------------------------\n" -" grt static syn -G -S[u] \n" +" grt static syn -G -S[u] -O \n" " [-M/[/]]\n" " [-T/////]\n" " [-F//] \n" +" [-X//] [-Y//]\n" " [-N] [-e] [-s]\n" -" -O \n" "\n" "\n\n" "Options:\n" @@ -137,6 +159,18 @@ printf("\n" " North, East and Vertical(Downward) Forces.\n" " Notice they will be scaled by .\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" +"\n" +" -Y//\n" +" Set the equidistant points in the east direction.\n" +" : start coordinate (km).\n" +" : end coordinate (km).\n" +" : sampling interval (km).\n" +"\n" " -N Components of results will be Z, N, E.\n" "\n" " -e Compute the spatial derivatives, ui_z and ui_r,\n" @@ -150,7 +184,7 @@ printf("\n" "Examples:\n" "----------------------------------------------------------------\n" " Say you have computed Static Green's functions with following command:\n" -" grt static greenfn -Mmilrow -D2/0 -X-5/5/10 -Y-5/5/10 -Ostgrn.nc\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" @@ -166,6 +200,9 @@ printf("\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" +"\n" +" You can also set a new XY grid, for example\n" +" grt static syn -Gstgrn.nc -Su1e16 -Ostsyn_ex.nc -X-5/5/0.5 -Y-5/5/0.5\n" "\n\n\n" "\n" ); @@ -179,7 +216,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:Nesh")) != -1) { + while ((opt = getopt(argc, argv, ":G:O:S:M:F:T:X:Y:Nesh")) != -1) { switch (opt) { // 输入 nc 文件名 case 'G': @@ -277,6 +314,52 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ } break; + // X坐标数组，-Xx1/x2/dx + case 'X': + Ctrl->X.active = true; + { + real_t a1, a2, delta; + if(3 != sscanf(optarg, "%lf/%lf/%lf", &a1, &a2, &delta)){ + GRTBadOptionError(X, ""); + }; + if(delta <= 0){ + GRTBadOptionError(X, "Can't set nonpositive dx(%f)", delta); + } + if(a1 > a2){ + GRTBadOptionError(X, "x1(%f) > x2(%f).", a1, a2); + } + + Ctrl->X.nx = floor((a2-a1)/delta) + 1; + Ctrl->X.xs = (real_t*)calloc(Ctrl->X.nx, sizeof(real_t)); + for(size_t i=0; iX.nx; ++i){ + Ctrl->X.xs[i] = a1 + delta*i; + } + } + break; + + // Y坐标数组，-Yy1/y2/dy + case 'Y': + Ctrl->Y.active = true; + { + real_t a1, a2, delta; + if(3 != sscanf(optarg, "%lf/%lf/%lf", &a1, &a2, &delta)){ + GRTBadOptionError(Y, ""); + }; + if(delta <= 0){ + GRTBadOptionError(Y, "Can't set nonpositive dy(%f)", delta); + } + if(a1 > a2){ + GRTBadOptionError(Y, "y1(%f) > y2(%f).", a1, a2); + } + + Ctrl->Y.ny = floor((a2-a1)/delta) + 1; + Ctrl->Y.ys = (real_t*)calloc(Ctrl->Y.ny, sizeof(real_t)); + for(size_t i=0; iY.ny; ++i){ + Ctrl->Y.ys[i] = a1 + delta*i; + } + } + break; + // 是否计算位移空间导数, 影响 calcUTypes 变量 case 'e': Ctrl->e.active = true; @@ -306,9 +389,142 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){ if(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 + if(Ctrl->X.active ^ Ctrl->Y.active){ + GRTRaiseError("If you want to set a new grid, you need set \"-X\" and \"-Y\" both.\n"); + } + Ctrl->isnewXYgrid = Ctrl->X.active; } +/** 选取震中距最近的点的格林函数计算新网格点的静态位移 */ +void grt_static_syn_new_xy( + size_t nx0, real_t *xs0, size_t ny0, real_t *ys0, + size_t nx, real_t *xs, size_t ny, real_t *ys, + realChnlGrid *u, realChnlGrid *uiz, realChnlGrid *uir, + GRT_SYN_TYPE computeType, real_t M0, real_t VpVs_ratio, 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 = nx0 * ny0; + size_t nr = nx * ny; + + // 原网格震中距序列 + real_t *rs0 = (real_t *)calloc(nr0, sizeof(real_t)); + for(size_t ix = 0; ix < nx0; ++ix){ + for(size_t iy = 0; iy < ny0; ++iy){ + rs0[iy + ix*ny0] = GRT_MAX(hypot(xs0[ix], ys0[iy]), GRT_MIN_DISTANCE); + } + } + + // 统计选取的最近震中距的绝对距离差的均值与标准差 + real_t mean_distDiff = 0.0; + real_t std_distDiff = 0.0; + real_t min_distDiff = __DBL_MAX__; + real_t max_distDiff = 0.0; + + realChnlGrid srcRadi = {0}; + + // 每个点逐个处理 + for(size_t ix = 0; ix < nx; ++ix){ + real_t x = xs[ix]; + for(size_t iy = 0; iy < ny; ++iy){ + real_t y = ys[iy]; + + // 震中距 + real_t dist = GRT_MAX(hypot(x, y), GRT_MIN_DISTANCE); + + // 方位角 + real_t azrad = atan2(y, x); + + // 从原震中距序列中找到最接近的 + size_t ir_pick = grt_findClosest_real_t(rs0, nr0, dist); + real_t distDiff = fabs(rs0[ir_pick] - dist); + mean_distDiff += distDiff; + std_distDiff += GRT_SQUARE(distDiff); + min_distDiff = GRT_MIN(min_distDiff, distDiff); + max_distDiff = GRT_MAX(max_distDiff, distDiff); + + size_t ir = iy + ix * ny; + + // 计算和位移相关量的种类（1-位移，2-ui_z，3-ui_r，4-ui_t） + int calcUTypes = (calc_upar)? 4 : 1; + real_t upar_scale = 1.0; + + realChnlGrid *up; // 使用对应类型的格林函数 + real_t tmpsyn[GRT_CHANNEL_NUM]; + + for(int ityp=0; ityp 0){ + switch (GRT_ZRT_CODES[ityp-1]){ + // 合成 ui_z, uir + case 'Z': case 'R': upar_scale = 1e-5; break; + // 合成 ui_t + case 'T': upar_scale = 1e-5 / dist; break; + default: break; + } + } + + if(ityp==1){ + up = uiz; + } else if(ityp==2){ + up = uir; + } else { + up = u; + } + + memset(tmpsyn, 0, sizeof(real_t)*GRT_CHANNEL_NUM); + + // 计算震源辐射因子 + grt_set_source_radiation(srcRadi, computeType, (ityp > 0) && GRT_ZRT_CODES[ityp-1]=='T', M0, upar_scale, VpVs_ratio, azrad, mchn); + + // 合成 + GRT_LOOP_ChnlGrid(im, c){ + int modr = GRT_SRC_M_ORDERS[im]; + if(modr==0 && GRT_ZRT_CODES[c]=='T') continue; + tmpsyn[c] += up[ir_pick][im][c] * srcRadi[im][c]; + } + + // 记录数据 + for(int i=0; iG.s_ingrid); @@ -383,11 +603,18 @@ int static_syn_main(int argc, char **argv){ NC_CHECK(nc_put_att_text(out_ncid, NC_GLOBAL, "computeType", strlen(Ctrl->s_computeType), Ctrl->s_computeType)); // 读入坐标变量 dimid, varid - size_t nx, ny; NC_CHECK(nc_inq_dimid(in_ncid, "north", &in_x_dimid)); - NC_CHECK(nc_inq_dimlen(in_ncid, in_x_dimid, &nx)); + NC_CHECK(nc_inq_dimlen(in_ncid, in_x_dimid, &nx0)); NC_CHECK(nc_inq_dimid(in_ncid, "east", &in_y_dimid)); - NC_CHECK(nc_inq_dimlen(in_ncid, in_y_dimid, &ny)); + NC_CHECK(nc_inq_dimlen(in_ncid, in_y_dimid, &ny0)); + xs0 = (real_t *)calloc(nx0, sizeof(real_t)); + ys0 = (real_t *)calloc(ny0, sizeof(real_t)); + + // 根据情况使用新网格 + nx = (Ctrl->isnewXYgrid)? Ctrl->X.nx : nx0; + ny = (Ctrl->isnewXYgrid)? Ctrl->Y.ny : ny0; + xs = (Ctrl->isnewXYgrid)? Ctrl->X.xs : xs0; + ys = (Ctrl->isnewXYgrid)? Ctrl->Y.ys : ys0; // 写入坐标变量 dimid, varid NC_CHECK(nc_def_dim(out_ncid, "north", nx, &out_x_dimid)); @@ -435,125 +662,64 @@ int static_syn_main(int argc, char **argv){ NC_CHECK(nc_enddef(out_ncid)); // 读取坐标变量 - real_t *xs = (real_t *)calloc(nx, sizeof(real_t)); - real_t *ys = (real_t *)calloc(ny, sizeof(real_t)); NC_CHECK(nc_inq_varid(in_ncid, "north", &in_x_varid)); - NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_x_varid, xs)); + NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_x_varid, xs0)); NC_CHECK(nc_inq_varid(in_ncid, "east", &in_y_varid)); - NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_y_varid, ys)); + NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_y_varid, ys0)); // 写入坐标变量 NC_CHECK(NC_FUNC_REAL(nc_put_var) (out_ncid, out_x_varid, xs)); NC_CHECK(NC_FUNC_REAL(nc_put_var) (out_ncid, out_y_varid, ys)); // 总震中距数 - size_t nr = nx * ny; + nr0 = nx0 * ny0; // 先将所有格林函数及其偏导读入内存， // 否则连续使用 nc_grt_var1 式读入效率太慢 - prealChnlGrid u; - prealChnlGrid uiz; - prealChnlGrid uir; - GRT_LOOP_ChnlGrid(im, c){ - int modr = GRT_SRC_M_ORDERS[im]; - // 先申请全 0 内存 - u[im][c] = (real_t *)calloc(nr, sizeof(real_t)); - uiz[im][c] = (real_t *)calloc(nr, sizeof(real_t)); - uir[im][c] = (real_t *)calloc(nr, sizeof(real_t)); - - if(modr==0 && GRT_ZRT_CODES[c]=='T') continue; + 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_get_var) (in_ncid, in_u_varids[im][c], u)); + for(size_t ir = 0; ir < nr0; ++ir){ + grn[ir][im][c] = u[ir]; + } - NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_u_varids[im][c], u[im][c])); + 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]; + } + } - if(Ctrl->e.active){ - NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_uiz_varids[im][c], uiz[im][c])); - NC_CHECK(NC_FUNC_REAL(nc_get_var) (in_ncid, in_uir_varids[im][c], uir[im][c])); } + GRT_SAFE_FREE_PTR(u); } + + // 新网格总震中距数 + nr = nx * ny; // 最终计算的结果 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); - // 每个点逐个处理 - for(size_t ix=0; ix < nx; ++ix){ - real_t x = xs[ix]; - for(size_t iy=0; iy < ny; ++iy){ - real_t y = ys[iy]; - - size_t ir = iy + ix*ny; - - // 方位角 - real_t azrad = atan2(y, x); - - // 震中距 - real_t dist = GRT_MAX(sqrt(x*x + y*y), GRT_MIN_DISTANCE); - - // 计算和位移相关量的种类（1-位移，2-ui_z，3-ui_r，4-ui_t） - int calcUTypes = (Ctrl->e.active)? 4 : 1; - real_t upar_scale = 1.0; - - real_t *(*up)[GRT_CHANNEL_NUM]; // 使用对应类型的格林函数 - real_t tmpsyn[GRT_CHANNEL_NUM]; - - for(int ityp=0; ityp 0){ - switch (GRT_ZRT_CODES[ityp-1]){ - // 合成 ui_z, uir - case 'Z': case 'R': upar_scale = 1e-5; break; - // 合成 ui_t - case 'T': upar_scale = 1e-5 / dist; break; - default: break; - } - } - - - if(ityp==1){ - up = uiz; - } else if(ityp==2){ - up = uir; - } else { - up = u; - } - - memset(tmpsyn, 0, sizeof(real_t)*GRT_CHANNEL_NUM); - - // 计算震源辐射因子 - grt_set_source_radiation(Ctrl->srcRadi, Ctrl->computeType, (ityp > 0) && GRT_ZRT_CODES[ityp-1]=='T', Ctrl->S.M0, upar_scale, src_va/src_vb, azrad, Ctrl->mchn); - - // 合成 - GRT_LOOP_ChnlGrid(im, c){ - int modr = GRT_SRC_M_ORDERS[im]; - if(modr==0 && GRT_ZRT_CODES[c]=='T') continue; - tmpsyn[c] += up[im][c][ir] * Ctrl->srcRadi[im][c]; - } - - // 记录数据 - for(int i=0; ie.active){ - grt_rot_zrt2zxy_upar(azrad, syn[ir], syn_upar[ir], dist*1e5); - } else { - grt_rot_zxy2zrt_vec(-azrad, syn[ir]); - } - } - - } - } + grt_static_syn_new_xy( + nx0, xs0, ny0, ys0, + nx, xs, ny, ys, + 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)); @@ -585,11 +751,12 @@ int static_syn_main(int argc, char **argv){ } // 释放内存 - GRT_LOOP_ChnlGrid(im, c){ - GRT_SAFE_FREE_PTR(u[im][c]); - GRT_SAFE_FREE_PTR(uiz[im][c]); - GRT_SAFE_FREE_PTR(uir[im][c]); - } + GRT_SAFE_FREE_PTR(xs0); + GRT_SAFE_FREE_PTR(ys0); + + 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); diff --git a/pygrt/c_interfaces.py b/pygrt/c_interfaces.py index 2c1ffa9..a3cf2fa 100755 --- a/pygrt/c_interfaces.py +++ b/pygrt/c_interfaces.py @@ -55,6 +55,20 @@ c_char_p ] +C_grt_static_syn_new_xy = libgrt.grt_static_syn_new_xy +"""合成新网格下的静态位移场""" +C_grt_static_syn_new_xy.restype = None +C_grt_static_syn_new_xy.argtypes = [ + c_size_t, PREAL, c_size_t, PREAL, + c_size_t, PREAL, c_size_t, PREAL, + POINTER((REAL*CHANNEL_NUM)*SRC_M_NUM), + POINTER((REAL*CHANNEL_NUM)*SRC_M_NUM), + POINTER((REAL*CHANNEL_NUM)*SRC_M_NUM), + c_int, REAL, REAL, REAL*MECHANISM_NUM, + c_bool, c_bool, + POINTER(REAL*CHANNEL_NUM), POINTER((REAL*CHANNEL_NUM)*CHANNEL_NUM), +] + C_grt_set_num_threads = libgrt.grt_set_num_threads """设置多线程数""" diff --git a/pygrt/pymod.py b/pygrt/pymod.py index a2c50de..e9d5655 100755 --- a/pygrt/pymod.py +++ b/pygrt/pymod.py @@ -559,8 +559,9 @@ def compute_grn( def compute_static_grn( self, - xarr:Union[np.ndarray,List[float],float], - yarr:Union[np.ndarray,List[float],float], + xarr:Union[np.ndarray,List[float],float,None]=None, + yarr:Union[np.ndarray,List[float],float,None]=None, + distarr:Union[np.ndarray,List[float],float,None]=None, keps:float=-1.0, k0:float=5.0, Length:float=15.0, @@ -572,10 +573,14 @@ def compute_static_grn( statsfile:Union[str,None]=None): r""" - Call the C function to calculate the static Green's functions and return them in a dict + Call the C function to calculate the static Green's functions and return them in a dict. + There're two ways to define the "epicentral distances": + 1. set both "xarr" and "yarr" to define a XY grid in advance. + 2. simply set "distarr", which equal to "xarr=[0.0], yarr=distarr". :param xarr: coordinate array in the north direction (km), or a single float. :param yarr: coordinate array in the east direction (km), or a single float. + :param distarr: equal to "xarr=[0.0], yarr=distarr" :param keps: automatic convergence condition, see (Yao and Harkrider (1983) for more details. negative value denotes not use. :param k0: k0 used to define the upper limit :math:`\tilde{k_{max}}=(k_{0}*\pi/hs)^2`, hs=max(abs(depsrc-deprcv),1.0) @@ -617,6 +622,19 @@ def compute_static_grn( depsrc = self.depsrc deprcv = self.deprcv + if distarr is not None: + if isinstance(distarr, float) or isinstance(distarr, int): + distarr = np.array([distarr*1.0]) + distarr = np.array(distarr) + + xarr = np.array([0.0]) + yarr = distarr.copy() + if np.any(yarr < 0.0): + raise ValueError("distances can't be negative.") + + if xarr is None or yarr is None: + raise ValueError("you need to set xarr and yarr or distarr.") + if isinstance(xarr, float) or isinstance(xarr, int): xarr = np.array([xarr*1.0]) xarr = np.array(xarr) diff --git a/pygrt/utils.py b/pygrt/utils.py index 051329b..7c76642 100755 --- a/pygrt/utils.py +++ b/pygrt/utils.py @@ -152,105 +152,94 @@ def _gen_syn_from_gf(st:Stream, calc_upar:bool, compute_type:GRT_SYN_TYPE, M0:fl return stall -def _gen_syn_from_static_gf(grn:dict, calc_upar:bool, compute_type:GRT_SYN_TYPE, M0:float, ZNE=False, **kwargs): +def _gen_syn_from_static_gf(grnDct:dict, calc_upar:bool, compute_type:GRT_SYN_TYPE, M0:float, ZNE=False, **kwargs): r""" 一个发起函数，根据不同震源参数，从静态格林函数中合成理论静态场 - :param grn: 计算好的静态格林函数, 字典类型 + :param grnDct: 计算好的静态格林函数, 字典类型 :param calc_upar: 是否计算位移u的空间导数 :param compute_type: 计算震源类型 :param M0: 标量地震矩, 单位dyne*cm :param ZNE: 是否以ZNE分量输出? + :param kwargs: 其它各种参数，包括震源参数，新网格参数等 """ chs = ZRTchs - sacin_prefixes = ["", "z", "r", ""] # 输入通道名 - srcName = ["EX", "VF", "HF", "DD", "DS", "SS"] - allchs = list(grn.keys()) + + mechn = _set_source_mechanism(compute_type, **kwargs) # 为张裂计算 Vp/Vs - src_va = grn['_src_va'] - src_vb = grn['_src_vb'] + src_va = grnDct['_src_va'] + src_vb = grnDct['_src_vb'] VpVs_ratio = src_va / src_vb - calcUTypes = 4 if calc_upar else 1 + xarr0:np.ndarray = grnDct['_xarr'] + yarr0:np.ndarray = grnDct['_yarr'] + nx0 = len(xarr0) + ny0 = len(yarr0) + nr0 = nx0 * ny0 - xarr:np.ndarray = grn['_xarr'] - yarr:np.ndarray = grn['_yarr'] + if "xarr" in kwargs and "yarr" in kwargs: + xarr:np.ndarray = kwargs['xarr'] + yarr:np.ndarray = kwargs['yarr'] + else: + xarr:np.ndarray = grnDct['_xarr'] + yarr:np.ndarray = grnDct['_yarr'] + + nx = len(xarr) + ny = len(yarr) + nr = nx * ny + + # 格林函数字典转为三个数组 + pygrn = np.zeros((nr0, SRC_M_NUM, CHANNEL_NUM), dtype=NPCT_REAL_TYPE, order='C'); c_pygrn = npct.as_ctypes(pygrn) + pygrn_uiz = np.zeros((nr0, SRC_M_NUM, CHANNEL_NUM), dtype=NPCT_REAL_TYPE, order='C'); c_pygrn_uiz = npct.as_ctypes(pygrn_uiz) + pygrn_uir = np.zeros((nr0, SRC_M_NUM, CHANNEL_NUM), dtype=NPCT_REAL_TYPE, order='C'); c_pygrn_uir = npct.as_ctypes(pygrn_uir) + if not calc_upar: + c_pygrn_uiz = c_pygrn_uir = None + for isrc in range(SRC_M_NUM): + src_name = SRC_M_NAME_ABBR[isrc] + for ic, comp in enumerate(ZRTchs): + pygrn[:,isrc,ic] = grnDct[f'{src_name}{comp}'].ravel() + if calc_upar: + pygrn_uiz[:,isrc,ic] = grnDct[f'z{src_name}{comp}'].ravel() + pygrn_uir[:,isrc,ic] = grnDct[f'r{src_name}{comp}'].ravel() + + # 结果数组 + syn = np.zeros((nr, CHANNEL_NUM), dtype=NPCT_REAL_TYPE, order='C'); c_syn = npct.as_ctypes(syn) + syn_upar = np.zeros((nr, CHANNEL_NUM, CHANNEL_NUM), dtype=NPCT_REAL_TYPE, order='C'); c_syn_upar = npct.as_ctypes(syn_upar) + + C_grt_static_syn_new_xy( + nx0, npct.as_ctypes(xarr0), ny0, npct.as_ctypes(yarr0), + nx, npct.as_ctypes(xarr), ny, npct.as_ctypes(yarr), + c_pygrn, c_pygrn_uiz, c_pygrn_uir, + compute_type.value, M0, VpVs_ratio, npct.as_ctypes(mechn), + ZNE, calc_upar, + c_syn, c_syn_upar, + ) # 结果字典 resDct = {} # 基本数据拷贝 - for k in grn.keys(): + for k in grnDct.keys(): if k[0] != '_': continue - resDct[k] = deepcopy(grn[k]) + resDct[k] = deepcopy(grnDct[k]) - XX = np.zeros((calcUTypes, 3, len(xarr), len(yarr)), dtype='f8') - dblsyn = (c_double*3)() - dblupar = (c_double*9)() + resDct['_xarr'] = xarr + resDct['_yarr'] = yarr - for iy in range(len(yarr)): - for ix in range(len(xarr)): - # 震中距 - dist = max(np.sqrt(xarr[ix]**2 + yarr[iy]**2), 1e-5) - - # 方位角 - azrad = np.arctan2(yarr[iy], xarr[ix]) - - upar_scale:float = 1.0 - for ityp in range(calcUTypes): - if ityp > 0: - upar_scale = 1e-5 - if ityp == 3: - upar_scale /= dist - - srcRadi = _set_source_radi(ityp==3, upar_scale, compute_type, M0, azrad, VpVs_ratio=VpVs_ratio, **kwargs) - - inpref = sacin_prefixes[ityp] - - for c in range(CHANNEL_NUM): - ch = chs[c] - for k in range(SRC_M_NUM): - coef = srcRadi[k, c] - if coef==0.0: - continue - - # 读入数据 - channel = f'{inpref}{srcName[k]}{ch}' - if channel not in allchs: - raise ValueError(f"Failed, channel=\"{channel}\" not exists.") - - XX[ityp, c, ix, iy] += coef*grn[channel][ix, iy] - - if ZNE: - for i in range(3): - dblsyn[i] = XX[0, i, ix, iy] - if calc_upar: - for k in range(3): - dblupar[k + i*3] = XX[i+1, k, ix, iy] - if calc_upar: - C_grt_rot_zrt2zxy_upar(azrad, dblsyn, dblupar, dist*1e5) - else: - C_grt_rot_zxy2zrt_vec(-azrad, dblsyn) - - for i in range(3): - XX[0, i, ix, iy] = dblsyn[i] - if calc_upar: - for k in range(3): - XX[i+1, k, ix, iy] = dblupar[k + i*3] - - - # 将XX数组分到字典中 if ZNE: chs = ZNEchs - for ityp in range(calcUTypes): - c1 = '' if ityp==0 else chs[ityp-1].lower() - for c in range(3): - resDct[f"{c1}{chs[c]}"] = XX[ityp, c].copy() - + for i1 in range(CHANNEL_NUM): + c1 = chs[i1] + resDct[c1] = syn[:, i1].reshape((nx, ny)) + if calc_upar: + for i2 in range(CHANNEL_NUM): + c2 = chs[i2] + resDct[f'{c2.lower()}{c1}'] = syn_upar[:, i2, i1].reshape((nx, ny)) + return resDct @@ -335,11 +324,34 @@ def _data_zrt2zne(stall:Stream): return stres -def _set_source_radi( - par_theta:bool, coef:float, compute_type:GRT_SYN_TYPE, M0:float, azrad:float, +def _set_source_mechanism( + compute_type:GRT_SYN_TYPE, fZ=None, fN=None, fE=None, strike=None, dip=None, rake=None, MT=None, **kwargs): + r""" + 整理 C 函数需要的震源机制数组 + """ + + mchn = np.zeros((MECHANISM_NUM,), dtype=NPCT_REAL_TYPE) + if compute_type == GRT_SYN_TYPE.GRT_SYN_EX: + pass + elif compute_type == GRT_SYN_TYPE.GRT_SYN_SF: + mchn[:3] = [fN, fE, fZ] + elif compute_type == GRT_SYN_TYPE.GRT_SYN_DC: + mchn[:3] = [strike, dip, rake] + elif compute_type == GRT_SYN_TYPE.GRT_SYN_TS: + mchn[:2] = [strike, dip] + elif compute_type == GRT_SYN_TYPE.GRT_SYN_MT: + mchn[:] = MT[:] + else: + raise ValueError("Unsupported source type.") + + return mchn + + +def _set_source_radi( + par_theta:bool, coef:float, compute_type:GRT_SYN_TYPE, M0:float, azrad:float, **kwargs): r""" 使用 C 函数计算不同震源的方向因子矩阵 @@ -358,20 +370,8 @@ def _set_source_radi( VpVs_ratio = kwargs['VpVs_ratio'] if 'VpVs_ratio' in kwargs else 0.0 src_coef = np.zeros((SRC_M_NUM, CHANNEL_NUM), dtype=NPCT_REAL_TYPE) - mchn = np.zeros((MECHANISM_NUM,), dtype=NPCT_REAL_TYPE) - if compute_type == GRT_SYN_TYPE.GRT_SYN_EX: - pass - elif compute_type == GRT_SYN_TYPE.GRT_SYN_SF: - mchn[:3] = [fN, fE, fZ] - elif compute_type == GRT_SYN_TYPE.GRT_SYN_DC: - mchn[:3] = [strike, dip, rake] - elif compute_type == GRT_SYN_TYPE.GRT_SYN_TS: - mchn[:2] = [strike, dip] - elif compute_type == GRT_SYN_TYPE.GRT_SYN_MT: - mchn[:] = MT[:] - else: - raise ValueError("Unsupported source type.") - + mchn = _set_source_mechanism(compute_type, **kwargs) + C_grt_set_source_radiation( npct.as_ctypes(src_coef), compute_type.value, par_theta, M0, coef, VpVs_ratio, azrad, npct.as_ctypes(mchn)) @@ -379,7 +379,7 @@ def _set_source_radi( return src_coef -def gen_syn_from_gf_DC(st:Union[Stream,dict], M0:float, strike:float, dip:float, rake:float, az:float=-999, ZNE=False, calc_upar:bool=False): +def gen_syn_from_gf_DC(st:Union[Stream,dict], M0:float, strike:float, dip:float, rake:float, az:float=-999, ZNE=False, calc_upar:bool=False, **kwargs): ''' Shear source, the unit of angles is all degrees(°) @@ -391,6 +391,8 @@ def gen_syn_from_gf_DC(st:Union[Stream,dict], M0:float, strike:float, dip:float, :param az: azimuth, 0 <= az <= 360 (not used for static case) :param ZNE: whether output in 'ZNE'-coord, default is 'ZRT' :param calc_upar: whether calculate the spatial derivatives of displacements. + :param kwargs: For static rsults, you can set "xarr" and "yarr" to define a new XY grid, + the program will automatically find the nearest Green's functions for each node. :return: - **stream** - :class:`obspy.Stream` @@ -400,11 +402,11 @@ def gen_syn_from_gf_DC(st:Union[Stream,dict], M0:float, strike:float, dip:float, raise ValueError(f"WRONG azimuth ({az})") return _gen_syn_from_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_DC, M0, az, ZNE, strike=strike, dip=dip, rake=rake) elif isinstance(st, dict): - return _gen_syn_from_static_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_DC, M0, ZNE, strike=strike, dip=dip, rake=rake) + return _gen_syn_from_static_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_DC, M0, ZNE, strike=strike, dip=dip, rake=rake, **kwargs) else: raise NotImplementedError -def gen_syn_from_gf_TS(st:Union[Stream,dict], M0:float, strike:float, dip:float, az:float=-999, ZNE=False, calc_upar:bool=False): +def gen_syn_from_gf_TS(st:Union[Stream,dict], M0:float, strike:float, dip:float, az:float=-999, ZNE=False, calc_upar:bool=False, **kwargs): ''' Tension source, the unit of angles is all degrees(°) @@ -415,6 +417,8 @@ def gen_syn_from_gf_TS(st:Union[Stream,dict], M0:float, strike:float, dip:float, :param az: azimuth, 0 <= az <= 360 (not used for static case) :param ZNE: whether output in 'ZNE'-coord, default is 'ZRT' :param calc_upar: whether calculate the spatial derivatives of displacements. + :param kwargs: For static rsults, you can set "xarr" and "yarr" to define a new XY grid, + the program will automatically find the nearest Green's functions for each node. :return: - **stream** - :class:`obspy.Stream` @@ -424,12 +428,12 @@ def gen_syn_from_gf_TS(st:Union[Stream,dict], M0:float, strike:float, dip:float, raise ValueError(f"WRONG azimuth ({az})") return _gen_syn_from_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_TS, M0, az, ZNE, strike=strike, dip=dip) elif isinstance(st, dict): - return _gen_syn_from_static_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_TS, M0, ZNE, strike=strike, dip=dip) + return _gen_syn_from_static_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_TS, M0, ZNE, strike=strike, dip=dip, **kwargs) else: raise NotImplementedError -def gen_syn_from_gf_SF(st:Union[Stream,dict], S:float, fN:float, fE:float, fZ:float, az:float=-999, ZNE=False, calc_upar:bool=False): +def gen_syn_from_gf_SF(st:Union[Stream,dict], S:float, fN:float, fE:float, fZ:float, az:float=-999, ZNE=False, calc_upar:bool=False, **kwargs): ''' Single-force source (dyne) @@ -441,6 +445,8 @@ def gen_syn_from_gf_SF(st:Union[Stream,dict], S:float, fN:float, fE:float, fZ:fl :param az: azimuth, 0 <= az <= 360 (not used for static case) :param ZNE: whether output in 'ZNE'-coord, default is 'ZRT' :param calc_upar: whether calculate the spatial derivatives of displacements. + :param kwargs: For static rsults, you can set "xarr" and "yarr" to define a new XY grid, + the program will automatically find the nearest Green's functions for each node. :return: - **stream** - :class:`obspy.Stream` @@ -450,12 +456,12 @@ def gen_syn_from_gf_SF(st:Union[Stream,dict], S:float, fN:float, fE:float, fZ:fl raise ValueError(f"WRONG azimuth ({az})") return _gen_syn_from_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_SF, S, az, ZNE, fN=fN, fE=fE, fZ=fZ) elif isinstance(st, dict): - return _gen_syn_from_static_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_SF, S, ZNE, fN=fN, fE=fE, fZ=fZ) + return _gen_syn_from_static_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_SF, S, ZNE, fN=fN, fE=fE, fZ=fZ, **kwargs) else: raise NotImplementedError -def gen_syn_from_gf_EX(st:Union[Stream,dict], M0:float, az:float=-999, ZNE=False, calc_upar:bool=False): +def gen_syn_from_gf_EX(st:Union[Stream,dict], M0:float, az:float=-999, ZNE=False, calc_upar:bool=False, **kwargs): ''' Explosion @@ -464,6 +470,8 @@ def gen_syn_from_gf_EX(st:Union[Stream,dict], M0:float, az:float=-999, ZNE=False :param az: azimuth, 0 <= az <= 360 (not used for static case) :param ZNE: whether output in 'ZNE'-coord, default is 'ZRT' :param calc_upar: whether calculate the spatial derivatives of displacements. + :param kwargs: For static rsults, you can set "xarr" and "yarr" to define a new XY grid, + the program will automatically find the nearest Green's functions for each node. :return: - **stream** - :class:`obspy.Stream` @@ -473,12 +481,12 @@ def gen_syn_from_gf_EX(st:Union[Stream,dict], M0:float, az:float=-999, ZNE=False raise ValueError(f"WRONG azimuth ({az})") return _gen_syn_from_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_EX, M0, az, ZNE) elif isinstance(st, dict): - return _gen_syn_from_static_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_EX, M0, ZNE) + return _gen_syn_from_static_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_EX, M0, ZNE, **kwargs) else: raise NotImplementedError -def gen_syn_from_gf_MT(st:Union[Stream,dict], M0:float, MT:ArrayLike, az:float=-999, ZNE=False, calc_upar:bool=False): +def gen_syn_from_gf_MT(st:Union[Stream,dict], M0:float, MT:ArrayLike, az:float=-999, ZNE=False, calc_upar:bool=False, **kwargs): ''' Moment tensor @@ -488,6 +496,8 @@ def gen_syn_from_gf_MT(st:Union[Stream,dict], M0:float, MT:ArrayLike, az:float=- :param az: azimuth, 0 <= az <= 360 (not used for static case) :param ZNE: whether output in 'ZNE'-coord, default is 'ZRT' :param calc_upar: whether calculate the spatial derivatives of displacements. + :param kwargs: For static rsults, you can set "xarr" and "yarr" to define a new XY grid, + the program will automatically find the nearest Green's functions for each node. :return: - **stream** - :class:`obspy.Stream` @@ -497,7 +507,7 @@ def gen_syn_from_gf_MT(st:Union[Stream,dict], M0:float, MT:ArrayLike, az:float=- raise ValueError(f"WRONG azimuth ({az})") return _gen_syn_from_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_MT, M0, az, ZNE, MT=MT) elif isinstance(st, dict): - return _gen_syn_from_static_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_MT, M0, ZNE, MT=MT) + return _gen_syn_from_static_gf(st, calc_upar, GRT_SYN_TYPE.GRT_SYN_MT, M0, ZNE, MT=MT, **kwargs) else: raise NotImplementedError diff --git a/test/_compare_c_py/compare.py b/test/_compare_c_py/compare.py index a2bd130..af74b84 100644 --- a/test/_compare_c_py/compare.py +++ b/test/_compare_c_py/compare.py @@ -1,129 +1,9 @@ import numpy as np import pygrt from obspy import * -from scipy.io import netcdf_file -from typing import List -import matplotlib.pyplot as plt +from compare_func import compare3, update_dict, static_compare3 -def compare3(st_py:Stream, c_prefix:str, ZNE:bool=False, dim2:bool=False): - """return average relative error""" - - if ZNE: - pattern = "[ZNE]" - else: - pattern = "[ZRT]" - - if dim2: - pattern = pattern*2 - - st_c = read(f"{c_prefix}{pattern}.sac") - - print(f"{c_prefix}{pattern}") - - error = 0.0 - nerr = 0 - for tr_c in st_c: - tr_py = st_py.select(channel=tr_c.stats.channel)[0] - if np.all(tr_c.data == 0.0) and np.all(tr_py.data == 0.0): - continue - - rerr = np.sum(np.abs(tr_c.data - tr_py.data)) / np.mean(np.abs(tr_py.data)) - print(tr_c.stats.channel, rerr) - - error += rerr - nerr += 1 - - - # fig, axs = plt.subplots(len(st_c), 1, figsize=(10, 1.5*len(st_c))) - # for i, tr_c in enumerate(st_c): - # ax = axs[i] - # t = np.arange(0, tr_c.stats.npts)*tr_c.stats.delta - # tr_py = st_py.select(channel=tr_c.stats.channel)[0] - # ax.plot(t, tr_c.data, label="c_"+tr_c.stats.channel) - # ax.plot(t, tr_py.data, label="py_"+tr_py.stats.channel) - # # 总误差水平线 - # ax.hlines(np.sum(np.abs(tr_c.data - tr_py.data)), *ax.get_xlim()) - # ax.legend() - # ax.grid() - - # plt.show() - - return error/nerr - -def update_dict(resDct:dict, Dct0:dict, prefix:str): - keys = resDct.keys() - for k in Dct0.keys(): - if k in keys: - continue - - resDct.update({f"{prefix}{k}": Dct0[k]}) - -def static_compare3(resDct:dict, c_prefix:str): - """return average relative error""" - - print(c_prefix) - - # read nc - with netcdf_file(c_prefix, mmap=False) as f: - error = 0.0 - nerr = 0 - - for k in f.variables: - if k == 'north' or k == 'east': - continue - val1 = resDct[k] - val2 = f.variables[k][:] - - if np.all(val1 == 0.0) and np.all(val2 == 0.0): - continue - - rerr = np.sum(np.abs(val1 - val2)) / np.mean(np.abs(val2)) - print(k, rerr) - - error += rerr - nerr += 1 - - # plot_static(resDct, c_prefix) - - return error/nerr - -def plot_static(resDct:dict, c_prefix:str): - n = len([k for k in resDct.keys() if k[0] != '_']) - fig, axs = plt.subplots(n, 3, figsize=(8, 3*n)) - with netcdf_file(c_prefix, mmap=False) as f: - norths = resDct['_xarr'] - easts = resDct['_yarr'] - - keys = f.variables - keys.pop('north') - keys.pop('east') - - for i, k in enumerate(keys): - val1 = resDct[k] - val2 = f.variables[k][:] - vmin = np.min(val1) - vmax = np.max(val1) - - norm = np.abs(val2) - norm[norm == 0.0] = 1.0 - - print(k, np.mean(np.abs(val1 - val2) / norm), np.max(np.abs(val1 - val2) / norm)) - - pcm0 = axs[i, 0].pcolorfast(easts, norths, val1.T, vmin=vmin, vmax=vmax) - pcm1 = axs[i, 1].pcolorfast(easts, norths, val2.T, vmin=vmin, vmax=vmax) - pcm2 = axs[i, 2].pcolorfast(easts, norths, ((val1 - val2) / norm).T) - - axs[i, 0].set_title("py_" + k) - axs[i, 1].set_title("c_" + k) - - fig.colorbar(pcm0) - fig.colorbar(pcm1) - fig.colorbar(pcm2) - - import os - fig.savefig(os.path.basename(c_prefix)[:-3] + ".pdf", bbox_inches='tight') - dist=10 depsrc=2 deprcv=3.3 diff --git a/test/_compare_c_py/compare_func.py b/test/_compare_c_py/compare_func.py new file mode 100644 index 0000000..fbef374 --- /dev/null +++ b/test/_compare_c_py/compare_func.py @@ -0,0 +1,124 @@ +import numpy as np +from obspy import * +from scipy.io import netcdf_file +from typing import List +import matplotlib.pyplot as plt + + +def compare3(st_py:Stream, c_prefix:str, ZNE:bool=False, dim2:bool=False): + """return average relative error""" + + if ZNE: + pattern = "[ZNE]" + else: + pattern = "[ZRT]" + + if dim2: + pattern = pattern*2 + + st_c = read(f"{c_prefix}{pattern}.sac") + + print(f"{c_prefix}{pattern}") + + error = 0.0 + nerr = 0 + for tr_c in st_c: + tr_py = st_py.select(channel=tr_c.stats.channel)[0] + if np.all(tr_c.data == 0.0) and np.all(tr_py.data == 0.0): + continue + + rerr = np.sum(np.abs(tr_c.data - tr_py.data)) / np.mean(np.abs(tr_py.data)) + print(tr_c.stats.channel, rerr) + + error += rerr + nerr += 1 + + + # fig, axs = plt.subplots(len(st_c), 1, figsize=(10, 1.5*len(st_c))) + # for i, tr_c in enumerate(st_c): + # ax = axs[i] + # t = np.arange(0, tr_c.stats.npts)*tr_c.stats.delta + # tr_py = st_py.select(channel=tr_c.stats.channel)[0] + # ax.plot(t, tr_c.data, label="c_"+tr_c.stats.channel) + # ax.plot(t, tr_py.data, label="py_"+tr_py.stats.channel) + # # 总误差水平线 + # ax.hlines(np.sum(np.abs(tr_c.data - tr_py.data)), *ax.get_xlim()) + # ax.legend() + # ax.grid() + + # plt.show() + + return error/nerr + +def update_dict(resDct:dict, Dct0:dict, prefix:str): + keys = resDct.keys() + for k in Dct0.keys(): + if k in keys: + continue + + resDct.update({f"{prefix}{k}": Dct0[k]}) + +def static_compare3(resDct:dict, c_prefix:str): + """return average relative error""" + + print(c_prefix) + + # read nc + with netcdf_file(c_prefix, mmap=False) as f: + error = 0.0 + nerr = 0 + + for k in f.variables: + if k == 'north' or k == 'east': + continue + val1 = resDct[k] + val2 = f.variables[k][:] + + if np.all(val1 == 0.0) and np.all(val2 == 0.0): + continue + + rerr = np.sum(np.abs(val1 - val2)) / np.mean(np.abs(val2)) + print(k, rerr) + + error += rerr + nerr += 1 + + # plot_static(resDct, c_prefix) + + return error/nerr + +def plot_static(resDct:dict, c_prefix:str): + n = len([k for k in resDct.keys() if k[0] != '_']) + fig, axs = plt.subplots(n, 3, figsize=(8, 3*n)) + with netcdf_file(c_prefix, mmap=False) as f: + norths = resDct['_xarr'] + easts = resDct['_yarr'] + + keys = f.variables + keys.pop('north') + keys.pop('east') + + for i, k in enumerate(keys): + val1 = resDct[k] + val2 = f.variables[k][:] + vmin = np.min(val1) + vmax = np.max(val1) + + norm = np.abs(val2) + norm[norm == 0.0] = 1.0 + + print(k, np.mean(np.abs(val1 - val2) / norm), np.max(np.abs(val1 - val2) / norm)) + + pcm0 = axs[i, 0].pcolorfast(easts, norths, val1.T, vmin=vmin, vmax=vmax) + pcm1 = axs[i, 1].pcolorfast(easts, norths, val2.T, vmin=vmin, vmax=vmax) + pcm2 = axs[i, 2].pcolorfast(easts, norths, ((val1 - val2) / norm).T) + + axs[i, 0].set_title("py_" + k) + axs[i, 1].set_title("c_" + k) + + fig.colorbar(pcm0) + fig.colorbar(pcm1) + fig.colorbar(pcm2) + + import os + fig.savefig(os.path.basename(c_prefix)[:-3] + ".pdf", bbox_inches='tight') diff --git a/test/_compare_c_py/compare_staticXY.py b/test/_compare_c_py/compare_staticXY.py new file mode 100644 index 0000000..9904d54 --- /dev/null +++ b/test/_compare_c_py/compare_staticXY.py @@ -0,0 +1,94 @@ +import numpy as np +import pygrt +from obspy import * +from compare_func import update_dict, static_compare3 + +depsrc=2 +deprcv=3.3 + +modname="../milrow" + +modarr = np.loadtxt(modname) + +pymod = pygrt.PyModel1D(modarr, depsrc, deprcv) + +#-------------------------- Static ----------------------------------------- +# 为了方便测试，避免引入其他因素的误差，这里有意避开 0 +xarr = np.arange(-3.1, 3.2, 0.6) +yarr = np.arange(-4.1, 4.2, 0.8) +S=1e24 + +fn=2 +fe=-1 +fz=4 + +stk=77 +dip=88 +rak=99 + +M11=1 +M12=-2 +M13=-5 +M22=0.5 +M23=3 +M33=1.2 + +static_grn = pymod.compute_static_grn(distarr=np.arange(0, 10+1e-8, 0.1), calc_upar=True) +AVGRERR2 = [] +# plot_static(static_grn, "static/stgrn.nc") + +for ZNE in [False, True]: + suffix = "-N" if ZNE else "" + static_syn = pygrt.utils.gen_syn_from_gf_EX(static_grn, S, ZNE=ZNE, calc_upar=True, xarr=xarr, yarr=yarr) + ststrain = pygrt.utils.compute_strain(static_syn) + strotation = pygrt.utils.compute_rotation(static_syn) + ststress = pygrt.utils.compute_stress(static_syn) + update_dict(static_syn, ststrain, "strain_") + update_dict(static_syn, ststress, "stress_") + update_dict(static_syn, strotation, "rotation_") + AVGRERR2.append(static_compare3(static_syn, f"static/stsyn_ex{suffix}.nc")) + + static_syn = pygrt.utils.gen_syn_from_gf_SF(static_grn, S, fn, fe, fz, ZNE=ZNE, calc_upar=True, xarr=xarr, yarr=yarr) + ststrain = pygrt.utils.compute_strain(static_syn) + strotation = pygrt.utils.compute_rotation(static_syn) + ststress = pygrt.utils.compute_stress(static_syn) + update_dict(static_syn, ststrain, "strain_") + update_dict(static_syn, ststress, "stress_") + update_dict(static_syn, strotation, "rotation_") + AVGRERR2.append(static_compare3(static_syn, f"static/stsyn_sf{suffix}.nc")) + + static_syn = pygrt.utils.gen_syn_from_gf_DC(static_grn, S, stk, dip, rak, ZNE=ZNE, calc_upar=True, xarr=xarr, yarr=yarr) + ststrain = pygrt.utils.compute_strain(static_syn) + strotation = pygrt.utils.compute_rotation(static_syn) + ststress = pygrt.utils.compute_stress(static_syn) + update_dict(static_syn, ststrain, "strain_") + update_dict(static_syn, ststress, "stress_") + update_dict(static_syn, strotation, "rotation_") + AVGRERR2.append(static_compare3(static_syn, f"static/stsyn_dc{suffix}.nc")) + + static_syn = pygrt.utils.gen_syn_from_gf_TS(static_grn, S, stk, dip, ZNE=ZNE, calc_upar=True, xarr=xarr, yarr=yarr) + ststrain = pygrt.utils.compute_strain(static_syn) + strotation = pygrt.utils.compute_rotation(static_syn) + ststress = pygrt.utils.compute_stress(static_syn) + update_dict(static_syn, ststrain, "strain_") + update_dict(static_syn, ststress, "stress_") + update_dict(static_syn, strotation, "rotation_") + AVGRERR2.append(static_compare3(static_syn, f"static/stsyn_ts{suffix}.nc")) + + static_syn = pygrt.utils.gen_syn_from_gf_MT(static_grn, S, [M11,M12,M13,M22,M23,M33], ZNE=ZNE, calc_upar=True, xarr=xarr, yarr=yarr) + ststrain = pygrt.utils.compute_strain(static_syn) + strotation = pygrt.utils.compute_rotation(static_syn) + ststress = pygrt.utils.compute_stress(static_syn) + update_dict(static_syn, ststrain, "strain_") + update_dict(static_syn, ststress, "stress_") + update_dict(static_syn, strotation, "rotation_") + AVGRERR2.append(static_compare3(static_syn, f"static/stsyn_mt{suffix}.nc")) + + +print("---------------- static --------------------") +AVGRERR2 = np.array(AVGRERR2) +print(AVGRERR2) +print(np.mean(AVGRERR2), np.min(AVGRERR2), np.max(AVGRERR2)) + +if np.mean(AVGRERR2) > 1e-5: + raise ValueError diff --git a/test/_compare_c_py/test_compare_staticXY_c_py.sh b/test/_compare_c_py/test_compare_staticXY_c_py.sh new file mode 100644 index 0000000..8bb9f1b --- /dev/null +++ b/test/_compare_c_py/test_compare_staticXY_c_py.sh @@ -0,0 +1,82 @@ +#!/bin/bash + +# 再测试静态解合成时传入新网格的结果 + +set -euo pipefail + +depsrc=2 +deprcv=3.3 + +modname="milrow" +out="GRN" + +rm -rf $out syn_* + + +#-------------------------- Static ----------------------------------------- +x1=-3.1 +x2=3.1 +dx=0.6 + +y1=-4.1 +y2=4.1 +dy=0.8 + +S=1e24 +fn=2 +fe=-1 +fz=4 + +stk=77 +dip=88 +rak=99 + +M11=1 +M12=-2 +M13=-5 +M22=0.5 +M23=3 +M33=1.2 + +rm -rf static + +mkdir -p static +cd static +grt static greenfn -M../../${modname} -D${depsrc}/${deprcv} -R0/10/0.1 -e -Ostgrn.nc + +for N in "" "-N" ; do +grt static syn -S$S -e $N -Gstgrn.nc -Ostsyn_ex$N.nc -X$x1/$x2/$dx -Y$y1/$y2/$dy +grt static strain stsyn_ex$N.nc +grt static rotation stsyn_ex$N.nc +grt static stress stsyn_ex$N.nc + +grt static syn -S$S -F$fn/$fe/$fz -e $N -Gstgrn.nc -Ostsyn_sf$N.nc -X$x1/$x2/$dx -Y$y1/$y2/$dy +grt static strain stsyn_sf$N.nc +grt static rotation stsyn_sf$N.nc +grt static stress stsyn_sf$N.nc + +grt static syn -S$S -M$stk/$dip/$rak -e $N -Gstgrn.nc -Ostsyn_dc$N.nc -X$x1/$x2/$dx -Y$y1/$y2/$dy +grt static strain stsyn_dc$N.nc +grt static rotation stsyn_dc$N.nc +grt static stress stsyn_dc$N.nc + +grt static syn -S$S -M$stk/$dip -e $N -Gstgrn.nc -Ostsyn_ts$N.nc -X$x1/$x2/$dx -Y$y1/$y2/$dy +grt static strain stsyn_ts$N.nc +grt static rotation stsyn_ts$N.nc +grt static stress stsyn_ts$N.nc + +grt static syn -S$S -T$M11/$M12/$M13/$M22/$M23/$M33 -e $N -Gstgrn.nc -Ostsyn_mt$N.nc -X$x1/$x2/$dx -Y$y1/$y2/$dy +grt static strain stsyn_mt$N.nc +grt static rotation stsyn_mt$N.nc +grt static stress stsyn_mt$N.nc +done + +cd - + + +# run pygrt and compare +python -u compare_staticXY.py + + +rm -rf $out syn_* +rm -rf static \ No newline at end of file diff --git a/test/_correct_static_disp/correctness.sh b/test/_correct_static_disp/correctness.sh index 1ec3c7a..96592f7 100755 --- a/test/_correct_static_disp/correctness.sh +++ b/test/_correct_static_disp/correctness.sh @@ -6,6 +6,7 @@ set -euo pipefail rm -rf *.nc # static greenfn +grt static greenfn -M../milrow -D2/3 -R0/6/0.2 -Cp -e -Ostgrn.nc grt static greenfn -M../milrow -D2/3 -X-3.1/3.1/0.2 -Y-2.1/2.1/0.2 -Cp -e -Ostgrn.nc # syn @@ -13,6 +14,7 @@ grt static syn -S1e20 -e -Gstgrn.nc -Ostsyn_ex.nc grt static syn -S1e20 -e -F2/-1/4 -Gstgrn.nc -Ostsyn_sf.nc grt static syn -S1e20 -e -M77/88/111 -Gstgrn.nc -Ostsyn_dc.nc grt static syn -S1e20 -e -M77/88 -Gstgrn.nc -Ostsyn_ts.nc +grt static syn -S1e20 -e -T1/-2/-5/0.5/3/1.2 -X-3.1/3.1/0.3 -Y-2.1/2.1/0.3 -Gstgrn.nc -Ostsyn_mt.nc # new XY grid grt static syn -S1e20 -e -T1/-2/-5/0.5/3/1.2 -Gstgrn.nc -Ostsyn_mt.nc # rotate to ZNE