Skip to content

Commit f5ac38b

Browse files
Dengda98Dengda98
authored
FEAT: update usage of -K, -L, -G, remove -V , in (static_)greenfn(#97)
1 parent 160ee7f commit f5ac38b

File tree

8 files changed

+268
-169
lines changed

8 files changed

+268
-169
lines changed

‎docs/source/Advanced/kernel/run/run.sh‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,8 +5,8 @@ rm -rf GRN* syn* pygrtstats* *.png
55
# -----------------------------------------------------------------
66
# BEGIN GRN
77
# -S-1 表示输出所有频率点的核函数
8-
# -V0.1 显式给定参考速度（用于定义波数积分上限），避免使用PTAM
8+
# -K+v0.1 显式给定参考速度（用于定义波数积分上限），避免使用PTAM
99
# -L20 定义波数积分间隔dk
10-
grt greenfn -Mmod1 -D0.01/0 -N500/0.02/0.8 -OGRN -R1 -V0.1 -S-1 -L20
10+
grt greenfn -Mmod1 -D0.01/0 -N500/0.02/0.8 -OGRN -R1 -K+v0.1 -S-1 -L20
1111
# END GRN
1212
# -----------------------------------------------------------------

‎example/far_field/run.sh‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,3 @@
11
#!/bin/bash
22

3-
grt greenfn -Mmilrow -OGRN -N1400/1 -D10/0 -R1800 -La20/1e-2/0 -E100 # -S-1 -G0/1/1/0
3+
grt greenfn -Mmilrow -OGRN -N1400/1 -D10/0 -R1800 -La20/1e-2/0 -E100 # -S-1 -Gvh

‎example/kernel_freq_response/run.sh‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,3 @@
11
#!/bin/bash
22

3-
grt greenfn -Mmod1 -D0.01/0 -N500/0.02/0.8 -OGRN -R1 -V0.1 -S-1 -L20
3+
grt greenfn -Mmod1 -D0.01/0 -N500/0.02/0.8 -OGRN -R1 -K+v0.1 -S-1 -L20

‎example/lamb_problem/run1.sh‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,4 +13,4 @@ modname="halfspace"
1313
out="GRN"
1414

1515
# compute Green's Functions
16-
grt greenfn -M${modname} -O${out} -N${nt}/${dt} -D${depsrc}/${deprcv} -R${dist} -G0/1/1/0 # Vertical and Horizontal Force
16+
grt greenfn -M${modname} -O${out} -N${nt}/${dt} -D${depsrc}/${deprcv} -R${dist} -Gvh # Vertical and Horizontal Force

‎example/multi_traces/run1.sh‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,4 +16,4 @@ modname="mod1"
1616
out="GRN"
1717

1818
# compute Green's Functions
19-
grt greenfn -M${modname} -O${out} -N${nt}/${dt} -D${depsrc}/${deprcv} -R${distarr} -G0/1/0/0 # 只生成垂直力源的格林函数
19+
grt greenfn -M${modname} -O${out} -N${nt}/${dt} -D${depsrc}/${deprcv} -R${distarr} -Gv # 只生成垂直力源的格林函数

‎pygrt/C_extension/include/grt/common/const.h‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -159,6 +159,13 @@ extern const char GRT_ZRT_CODES[];
159159
extern const char GRT_ZNE_CODES[];
160160

161161

162+
/** 波数积分方法 */
163+
enum {
164+
GRT_K_INTEG_METHOD_DWM = 0, // 离散波数法
165+
GRT_K_INTEG_METHOD_FIM, // 固定间隔 Filon 积分
166+
GRT_K_INTEG_METHOD_SAFIM, // 自适应 Filon 积分
167+
};
168+
162169
/**
163170
* 设置OpenMP多线程数
164171
*

‎pygrt/C_extension/src/dynamic/grt_greenfn.c‎

Lines changed: 146 additions & 95 deletions
Large diffs are not rendered by default.

‎pygrt/C_extension/src/static/grt_static_greenfn.c‎

Lines changed: 109 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@
1919
#include "grt.h"
2020

2121
// 一些变量的非零默认值
22-
#define GRT_GREENFN_V_VMIN_REF 0.1
22+
#define GRT_GREENFN_K_VMIN 0.1
2323
#define GRT_GREENFN_K_K0 5.0
2424
#define GRT_GREENFN_L_LENGTH 15.0
2525

@@ -45,6 +45,7 @@ typedef struct {
4545
/** 波数积分间隔 */
4646
struct {
4747
bool active;
48+
MYINT method;
4849
MYREAL Length;
4950
MYREAL filonLength;
5051
MYREAL safilonTol;
@@ -55,12 +56,9 @@ typedef struct {
5556
bool active;
5657
MYREAL keps;
5758
MYREAL k0;
59+
MYREAL vmin;
60+
bool v_active;
5861
} K;
59-
/** 参考速度 */
60-
struct {
61-
bool active;
62-
MYREAL vmin_ref;
63-
} V;
6462
/** 波数积分过程的核函数文件 */
6563
struct {
6664
bool active;
@@ -131,8 +129,8 @@ printf("\n"
131129
"Usage:\n"
132130
"----------------------------------------------------------------\n"
133131
" grt static greenfn -M<model> -D<depsrc>/<deprcv> -X<x1>/<x2>/<dx> \n"
134-
" -Y<y1>/<y2>/<dy> [-L<length>] [-V<vmin_ref>] \n"
135-
" [-K<k0>[/<keps>]] [-S] [-e]\n"
132+
" -Y<y1>/<y2>/<dy> [-L<length>] \n"
133+
" [-K[+k<k0>][+e<keps>][+v<vmin>]] [-S] [-e]\n"
136134
"\n\n"
137135
"Options:\n"
138136
"----------------------------------------------------------------\n"
@@ -161,37 +159,24 @@ printf("\n"
161159
" <y2>: end coordinate (km).\n"
162160
" <dy>: sampling interval (km).\n"
163161
"\n"
164-
" -L[a]<length>[/<Flength>/<Fcut>]\n"
162+
" -L[a|l]<length>[/<Flength>|<Ftol>/<Fcut>]\n"
165163
" Define the wavenumber integration interval\n"
166164
" dk=(2*PI)/(<length>*rmax). rmax is the maximum \n"
167165
" epicentral distance. \n"
168-
" There are 3 cases:\n"
166+
" There are 4 cases:\n"
169167
" + (default) not set or set 0.0.\n"); printf(
170168
" <length> will be %.1f.\n", GRT_GREENFN_L_LENGTH); printf(
171-
" + manually set one POSITIVE value, e.g. -L20\n"
169+
" + manually set one POSITIVE value, e.g. -Ll20\n"
172170
" + manually set three POSITIVE values, \n"
173-
" e.g. -L20/5/10, means split the integration \n"
171+
" e.g. -Ll20/5/10, means split the integration \n"
174172
" into two parts, [0, k*] and [k*, kmax], \n"
175173
" in which k*=<Fcut>/rmax, and use DWM with\n"
176174
" <length> and FIM with <Flength>, respectively.\n"
177175
" + manually set three POSITIVE values, with -La,\n"
178176
" in this case, <Flength> will be <Ftol> for Self-\n"
179177
" Adaptive FIM.\n"
180178
"\n"
181-
" -V<vmin_ref> \n"
182-
" (Inherited from the dynamic case, and the numerical\n"
183-
" value will not be used in here, except its sign.)\n"
184-
" + (default) not set or set 0.0.\n"); printf(
185-
" <vmin_ref> will be the minimum velocity\n"
186-
" of model, but limited to %.1f. and if the \n", GRT_GREENFN_V_VMIN_REF); printf(
187-
" depth gap between source and receiver is \n"
188-
" thinner than %.1f km, PTAM will be appled\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf(
189-
" automatically.\n"
190-
" + manually set POSITIVE value. \n"
191-
" + manually set NEGATIVE value, \n"
192-
" and PTAM will be appled.\n"
193-
"\n"
194-
" -K<k0>[/<keps>]\n"
179+
" -K[+k<k0>][+e<keps>][+v<vmin>]\n"
195180
" Several parameters designed to define the\n"
196181
" behavior in wavenumber integration. The upper\n"
197182
" bound is k0,\n"
@@ -201,7 +186,19 @@ printf("\n"
201186
" <keps>: a threshold for break wavenumber \n"
202187
" integration in advance. See \n"
203188
" (Yao and Harkrider, 1983) for details.\n"
204-
" Default 0.0 not use.\n"); printf(
189+
" Default 0.0 not use.\n"
190+
" <vmin>: Minimum velocity (km/s) for reference. This\n"
191+
" is designed to define the upper bound \n"
192+
" of wavenumber integration.\n"
193+
" There are 3 cases:\n"
194+
" + (default) not set or set 0.0.\n"); printf(
195+
" <vmin> will be the minimum velocity\n"
196+
" of model, but limited to %.1f. and if \n", GRT_GREENFN_K_VMIN); printf(
197+
" hs is thinner than %.1f km, PTAM will be appled\n", GRT_MIN_DEPTH_GAP_SRC_RCV); printf(
198+
" automatically.\n"
199+
" + manually set POSITIVE value. \n"
200+
" + manually set NEGATIVE value, \n"
201+
" and PTAM will be appled.\n"
205202
"\n"
206203
" -S Output statsfile in wavenumber integration.\n"
207204
"\n"
@@ -232,7 +229,7 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
232229
char* command = Ctrl->name;
233230

234231
// 先为个别参数设置非0初始值
235-
Ctrl->V.vmin_ref = GRT_GREENFN_V_VMIN_REF;
232+
Ctrl->K.vmin = GRT_GREENFN_K_VMIN;
236233
Ctrl->K.k0 = GRT_GREENFN_K_K0;
237234

238235
int opt;
@@ -266,50 +263,94 @@ static void getopt_from_command(GRT_MODULE_CTRL *Ctrl, int argc, char **argv){
266263
}
267264
break;
268265

269-
// 波数积分间隔 -L[a]<length>[/<Flength>/<Fcut>]
266+
// 波数积分间隔 -L[a|l]<length>[/<Flength>|<Ftol>/<Fcut>]
270267
case 'L':
271268
Ctrl->L.active = true;
272-
{
273-
// 检查首字母是否为a，表明使用自适应Filon积分
274-
int pos=0;
275-
bool useSAFIM = false;
276-
if(optarg[0] == 'a'){
277-
pos++;
278-
useSAFIM = true;
269+
{
270+
// 若是纯数字，即未指定子选项，则直接使用DWM，并指定步长
271+
if(isdigit(optarg[0])){
272+
Ctrl->L.method = GRT_K_INTEG_METHOD_DWM;
273+
// 仅接受一个值，若有多个值的分隔符则报错
274+
if(strchr(optarg, '/') != NULL){
275+
GRTBadOptionError(command, L, "single -L accept only 1 argument, but found %s.\n", optarg);
276+
}
277+
if(1 != sscanf(optarg, "%lf", &Ctrl->L.Length)){
278+
GRTBadOptionError(command, L, "");
279+
}
279280
}
280-
double filona = 0.0;
281-
int n = sscanf(optarg+pos, "%lf/%lf/%lf", &Ctrl->L.Length, &filona, &Ctrl->L.filonCut);
282-
if(n != 1 && n != 3){
283-
GRTBadOptionError(command, L, "");
284-
};
285-
if(n == 1 && Ctrl->L.Length <= 0){
286-
GRTBadOptionError(command, L, "Length should be positive.");
287-
}
288-
if(n == 3 && (filona <= 0 || Ctrl->L.filonCut < 0)){
289-
GRTBadOptionError(command, L, "Flength/Ftol should be positive, Fcut should be nonnegative.");
290-
}
291-
if(n == 3){
292-
useSAFIM ? (Ctrl->L.safilonTol = filona) : (Ctrl->L.filonLength = filona);
281+
// 指定了子选项
282+
else {
283+
// 固定间隔 Filon 积分
284+
if(optarg[0] == 'l'){
285+
Ctrl->L.method = GRT_K_INTEG_METHOD_FIM;
286+
if(3 != sscanf(optarg+1, "%lf/%lf/%lf", &Ctrl->L.Length, &Ctrl->L.filonLength, &Ctrl->L.filonCut)){
287+
GRTBadOptionError(command, L, "");
288+
}
289+
if(Ctrl->L.filonLength <= 0.0){
290+
GRTBadOptionError(command, L, "Flength should be positive.");
291+
}
292+
}
293+
// 自适应采样
294+
else if(optarg[0] == 'a'){
295+
Ctrl->L.method = GRT_K_INTEG_METHOD_SAFIM;
296+
if(3 != sscanf(optarg+1, "%lf/%lf/%lf", &Ctrl->L.Length, &Ctrl->L.safilonTol, &Ctrl->L.filonCut)){
297+
GRTBadOptionError(command, L, "");
298+
}
299+
if(Ctrl->L.safilonTol <= 0.0){
300+
GRTBadOptionError(command, L, "safilonTol should be positive.");
301+
}
302+
}
303+
304+
// 检查共有参数
305+
if(Ctrl->L.Length <= 0.0){
306+
GRTBadOptionError(command, L, "Length should be positive.");
307+
}
308+
if(Ctrl->L.filonCut < 0.0){
309+
GRTBadOptionError(command, L, "Fcut should be nonnegative.");
310+
}
293311
}
294312
}
295313
break;
296314

297-
// 参考最小速度 -Vvmin_ref
298-
case 'V':
299-
Ctrl->V.active = true;
300-
if(0 == sscanf(optarg, "%lf", &Ctrl->V.vmin_ref)){
301-
GRTBadOptionError(command, V, "");
302-
};
303-
break;
304-
305-
// 波数积分相关变量 -Kk0/keps
315+
// 波数积分相关变量 -K[+k<k0>][+e<keps>][+v<vmin>]
306316
case 'K':
307317
Ctrl->K.active = true;
308-
if(0 == sscanf(optarg, "%lf/%lf", &Ctrl->K.k0, &Ctrl->K.keps)){
309-
GRTBadOptionError(command, K, "");
310-
};
311-
if(Ctrl->K.k0 < 0.0){
312-
GRTBadOptionError(command, K, "Can't set negative k0(%f).", Ctrl->K.k0);
318+
{
319+
char *line = strdup(optarg);
320+
char *token = strtok(line, "+");
321+
while(token != NULL){
322+
switch(token[0]) {
323+
case 'k':
324+
if(1 != sscanf(token+1, "%lf", &Ctrl->K.k0)){
325+
GRTBadOptionError(command, K+k, "");
326+
}
327+
if(Ctrl->K.k0 < 0.0){
328+
GRTBadOptionError(command, K, "Can't set negative k0(%f).", Ctrl->K.k0);
329+
}
330+
break;
331+
332+
case 'e':
333+
if(1 != sscanf(token+1, "%lf", &Ctrl->K.keps)){
334+
GRTBadOptionError(command, K+e, "");
335+
}
336+
break;
337+
338+
case 'v':
339+
Ctrl->K.v_active = true;
340+
if(1 != sscanf(token+1, "%lf", &Ctrl->K.vmin)){
341+
GRTBadOptionError(command, K+v, "");
342+
}
343+
break;
344+
345+
default:
346+
GRTBadOptionError(command, K, "-K+%s is not supported.", token);
347+
break;
348+
}
349+
350+
token = strtok(NULL, "+");
351+
}
352+
353+
GRT_SAFE_FREE_PTR(line);
313354
}
314355
break;
315356

@@ -453,13 +494,13 @@ int static_greenfn_main(int argc, char **argv){
453494
grt_get_mod1d_vmin_vmax(mod1d, &vmin, &vmax);
454495

455496
// 参考最小速度
456-
if(!Ctrl->V.active){
457-
Ctrl->V.vmin_ref = GRT_MAX(vmin, GRT_GREENFN_V_VMIN_REF);
497+
if(!Ctrl->K.v_active){
498+
Ctrl->K.vmin = GRT_MAX(vmin, GRT_GREENFN_K_VMIN);
458499
}
459500

460501
// 如果没有主动设置vmin_ref，则判断是否要自动使用PTAM
461-
if( !Ctrl->V.active && fabs(Ctrl->D.deprcv - Ctrl->D.depsrc) <= GRT_MIN_DEPTH_GAP_SRC_RCV) {
462-
Ctrl->V.vmin_ref = - fabs(Ctrl->V.vmin_ref);
502+
if( !Ctrl->K.v_active && fabs(Ctrl->D.deprcv - Ctrl->D.depsrc) <= GRT_MIN_DEPTH_GAP_SRC_RCV) {
503+
Ctrl->K.vmin = - fabs(Ctrl->K.vmin);
463504
}
464505

465506
// 设置积分间隔默认值
@@ -484,7 +525,7 @@ int static_greenfn_main(int argc, char **argv){
484525
//==============================================================================
485526
// 计算静态格林函数
486527
grt_integ_static_grn(
487-
mod1d, Ctrl->nr, Ctrl->rs, Ctrl->V.vmin_ref, Ctrl->K.keps, Ctrl->K.k0, Ctrl->L.Length, Ctrl->L.filonLength, Ctrl->L.safilonTol, Ctrl->L.filonCut,
528+
mod1d, Ctrl->nr, Ctrl->rs, Ctrl->K.vmin, Ctrl->K.keps, Ctrl->K.k0, Ctrl->L.Length, Ctrl->L.filonLength, Ctrl->L.safilonTol, Ctrl->L.filonCut,
488529
grn, Ctrl->e.active, grn_uiz, grn_uir,
489530
Ctrl->S.s_statsdir
490531
);

0 commit comments

Comments
 (0)