Skip to content

Commit 152906f

Browse files
Dengda98Dengda98
committed
add slowness-shape check, add a way to compute traveltime seg-by-seg
1 parent 57da6ce commit 152906f

File tree

6 files changed

+59
-9
lines changed

6 files changed

+59
-9
lines changed

‎.gitignore‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,4 +8,5 @@
88
/pyfmm/__pycache__
99
/pyfmm/C_extension/build
1010
/pyfmm/C_extension/lib
11+
/pyfmm/C_extension/bin
1112
/test_w

‎pyfmm/C_extension/bin/grt.fmm‎

-209 KB
Binary file not shown.

‎pyfmm/C_extension/include/fmm.h‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -210,6 +210,7 @@ MYREAL get_neighbour_travt(
210210
* @param pp (in)接收点维度3坐标
211211
* @param seglen (in)射线段长度
212212
* @param segfac (in)t < segfac*seglen/v，当射线追踪到在源点附近时，射线直接连接源点
213+
* @param Slw (in)展平的三维慢度场，若非NULL则使用累加求和计算走时，否则直接从走时场中插值得到走时
213214
* @param TT (in)展平的三维走时场
214215
* @param sphcoord (in)是否使用球坐标
215216
* @param rays (out)输出展平的三维射线
@@ -224,5 +225,5 @@ MYREAL FMM_raytracing(
224225
const double *ps, int np,
225226
double r0, double t0, double p0,
226227
double rr, double tt, double pp, double seglen, double segfac,
227-
const MYREAL *TT, bool sphcoord,
228+
const MYREAL *Slw, const MYREAL *TT, bool sphcoord,
228229
double *rays, int *N);

‎pyfmm/C_extension/src/fmm.c‎

Lines changed: 28 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -810,14 +810,15 @@ MYREAL FMM_raytracing(
810810
const double *ps, int np,
811811
double r0, double t0, double p0,
812812
double rr, double tt, double pp, double seglen, double segfac,
813-
const MYREAL *TT, bool sphcoord,
813+
const MYREAL *Slw, const MYREAL *TT, bool sphcoord,
814814
// MYREAL *gTr, MYREAL *gTt, MYREAL *gTp,
815815
double *rays, int *N)
816816
{
817817
double dr = (nr>1)? rs[1] - rs[0] : 1e-6;
818818
double dt = (nt>1)? ts[1] - ts[0] : 1e-6;
819819
double dp = (np>1)? ps[1] - ps[0] : 1e-6;
820820
double seglen0 = seglen;
821+
double seglen1;
821822

822823
double xx, yy, zz, x0, y0, z0;
823824
double dx, dy, dz;
@@ -843,6 +844,7 @@ MYREAL FMM_raytracing(
843844
double r1, t1, p1;
844845
double x1, y1, z1;
845846
double r11, t11, p11;
847+
double rmid, tmid, pmid, vmid;
846848
r1 = rr;
847849
t1 = tt;
848850
p1 = pp;
@@ -871,6 +873,7 @@ MYREAL FMM_raytracing(
871873
rs, nr, ts, nt, ps, np, ntp, TT, r1, t1, p1,
872874
&gtr, &gtt, &gtp, NULL, NULL);
873875
MYREAL trem = travt, trem1;
876+
MYREAL travt1 = 0.0;
874877

875878
// normalize gradient
876879
gtr /= dr;
@@ -903,6 +906,7 @@ MYREAL FMM_raytracing(
903906
// update
904907
seglen = seglen0;
905908
for(int i=0; i<5; ++i){
909+
seglen1 = seglen;
906910
if(sphcoord){
907911
p11 = p1 - gtp*seglen/(r1*sin(t1));
908912
t11 = t1 - gtt*seglen/r1;
@@ -959,6 +963,17 @@ MYREAL FMM_raytracing(
959963
gtt /= norm;
960964
gtp /= norm;
961965

966+
// 求中点速度，算走时
967+
// 这里做了近似，认为射线段比较小，不论坐标系为何，直接去平均值作为中点
968+
if(Slw != NULL){
969+
rmid = (r1 + rays[3*idot-3])/2.0;
970+
tmid = (t1 + rays[3*idot-2])/2.0;
971+
pmid = (p1 + rays[3*idot-1])/2.0;
972+
travt1 += trilinear_one_ravel(
973+
rs, nr, ts, nt, ps, np, ntp, Slw, rmid, tmid, pmid,
974+
NULL, NULL, NULL, NULL, NULL) * seglen1;
975+
}
976+
962977
} // END tracing
963978

964979
rays[3*idot] = r0;
@@ -968,7 +983,18 @@ MYREAL FMM_raytracing(
968983
idot++;
969984
*N = idot;
970985

971-
return travt;
986+
if(Slw != NULL){
987+
rmid = (r0 + rays[3*idot-3])/2.0;
988+
tmid = (t0 + rays[3*idot-2])/2.0;
989+
pmid = (p0 + rays[3*idot-1])/2.0;
990+
travt1 += trilinear_one_ravel(
991+
rs, nr, ts, nt, ps, np, ntp, Slw, rmid, tmid, pmid,
992+
NULL, NULL, NULL, NULL, NULL) * seglen1;
993+
return travt1;
994+
} else {
995+
return travt;
996+
}
997+
972998
}
973999

9741000

‎pyfmm/c_interfaces.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -72,7 +72,7 @@ def load_c_lib(use_float:bool=False):
7272
PDOUBLE, c_int,
7373
c_double, c_double, c_double,
7474
c_double, c_double, c_double, c_double, c_double,
75-
PREAL, c_bool,
75+
PREAL, PREAL, c_bool,
7676
PDOUBLE, PINT
7777
]
7878

‎pyfmm/traveltime.py‎

Lines changed: 27 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,7 @@ def travel_time_source(
5656
global FSM_nsweep
5757

5858
check_xyz_arr(xarr, yarr, zarr, sphcoord)
59-
check_slowness(slw)
59+
check_slowness(xarr, yarr, zarr, slw)
6060

6161
# 对于并行情况，至少迭代两次
6262
if FSMparallel:
@@ -139,7 +139,7 @@ def travel_time_iniTT(
139139
global FSM_nsweep
140140

141141
check_xyz_arr(xarr, yarr, zarr, sphcoord)
142-
check_slowness(slw)
142+
check_slowness(xarr, yarr, zarr, slw)
143143

144144
# 对于并行情况，至少迭代两次
145145
if FSMparallel:
@@ -181,7 +181,7 @@ def travel_time_iniTT(
181181
def raytracing(
182182
TT:np.ndarray, srcloc:list, rcvloc:list,
183183
xarr:np.ndarray, yarr:np.ndarray, zarr:np.ndarray,
184-
seglen:float, segfac:int=3, sphcoord:bool=False, maxdots:int=10000):
184+
seglen:float, slw:np.ndarray|None=None, segfac:int=3, sphcoord:bool=False, maxdots:int=10000):
185185
r'''
186186
根据给定源点坐标计算的走时场，使用梯度下降法做射线追踪
187187
@@ -192,13 +192,17 @@ def raytracing(
192192
:param yarr: :math:`y` 或 :math:`\theta` 节点坐标数组，要求等距升序排列
193193
:param zarr: :math:`z` 或 :math:`\phi` 节点坐标数组，要求等距升序排列
194194
:param seglen: 射线段长度，与xyz的长度量纲保持一致
195+
:param slw: 形状为(nx, ny, nz)的三维慢度场，若非None则使用累加求和计算走时，否则直接从走时场中插值得到走时
195196
:param segfac: t < segfac*seglen/v，当射线追踪到在源点附近时，射线直接连接源点
196197
:param sphcoord: 是否使用球坐标
197198
:param maxdots: 射线最大点数
198199
199200
:return: (接收点走时，形状为(ndots, 3)的射线坐标)
200201
'''
201202

203+
if isinstance(slw, np.ndarray):
204+
check_slowness(xarr, yarr, zarr, slw)
205+
202206
sx, sy, sz = np.array(srcloc).astype('f8')
203207
rx, ry, rz = np.array(rcvloc).astype('f8')
204208

@@ -225,6 +229,10 @@ def raytracing(
225229

226230
TT_ravel = TT.ravel().astype(c_interfaces.NPCT_REAL_TYPE)
227231
c_TT = npct.as_ctypes(TT_ravel)
232+
c_slw = None
233+
if slw is not None:
234+
slw_ravel = slw.ravel().astype(c_interfaces.NPCT_REAL_TYPE)
235+
c_slw = npct.as_ctypes(slw_ravel)
228236

229237
c_xarr = npct.as_ctypes(xarr.astype('f8'))
230238
c_yarr = npct.as_ctypes(yarr.astype('f8'))
@@ -240,7 +248,7 @@ def raytracing(
240248
c_zarr, len(zarr),
241249
sx, sy, sz,
242250
rx, ry, rz, float(seglen), int(segfac),
243-
c_TT, sphcoord,
251+
c_slw, c_TT, sphcoord,
244252
c_rays, byref(c_ndots)
245253
)
246254

@@ -278,6 +286,14 @@ def check_xyz_arr(
278286
检查三个维度的数组是否符合基本要求
279287
'''
280288

289+
# 检查类型
290+
if not isinstance(xarr, np.ndarray):
291+
raise ValueError("xarr should be an instance of numpy.ndarray.")
292+
if not isinstance(yarr, np.ndarray):
293+
raise ValueError("yarr should be an instance of numpy.ndarray.")
294+
if not isinstance(zarr, np.ndarray):
295+
raise ValueError("zarr should be an instance of numpy.ndarray.")
296+
281297
# 检查维数
282298
if len(xarr)==0:
283299
raise ValueError("xarr is empty.")
@@ -308,11 +324,17 @@ def check_xyz_arr(
308324

309325

310326

311-
def check_slowness(slw:np.ndarray):
327+
def check_slowness(xarr:np.ndarray, yarr:np.ndarray, zarr:np.ndarray, slw:np.ndarray):
312328
r'''
313329
对慢度数组进行检查
314330
'''
331+
if not isinstance(slw, np.ndarray):
332+
raise ValueError("slw should be an instance of numpy.ndarray.")
333+
315334
if np.any(slw <= 0.0):
316335
raise ValueError("Slowness should be positive.")
317336

337+
shapexyz = (len(xarr), len(yarr), len(zarr))
338+
if slw.shape != shapexyz:
339+
raise ValueError(f"Shape of slowness should be {shapexyz}, but {slw.shape}.")
318340

0 commit comments

Comments
 (0)