Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 24 additions & 10 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
cmake_minimum_required(VERSION 3.18)

project( SurfATT VERSION 2.1.0 LANGUAGES CXX Fortran)
project( SurfATT VERSION 2.1.1 LANGUAGES CXX Fortran)

include(GNUInstallDirs)

Expand Down Expand Up @@ -143,9 +143,12 @@ endif()
add_subdirectory(external_libs)

# yaml-cpp: prefer a system installation, fall back to the bundled submodule.
find_package(yaml-cpp QUIET)
if (yaml-cpp_FOUND AND NOT ${USE_EXTERNAL_LIBS})
message(STATUS "Found yaml-cpp: ${YAML_CPP_INCLUDE_DIR}")
if (NOT USE_EXTERNAL_LIBS)
find_package(yaml-cpp QUIET)
endif()

if (yaml-cpp_FOUND AND NOT USE_EXTERNAL_LIBS)
message(STATUS "Found system yaml-cpp: ${yaml-cpp_DIR}")
# Normalize to yaml-cpp::yaml-cpp to avoid deprecated target warnings.
if (NOT TARGET yaml-cpp::yaml-cpp)
if (TARGET yaml-cpp)
Expand All @@ -162,17 +165,28 @@ if (yaml-cpp_FOUND AND NOT ${USE_EXTERNAL_LIBS})
endif()
endif()
endif()
else()
message(STATUS "System yaml-cpp not found (or bundled libraries requested); using external_libs/yaml-cpp")
initialize_submodule(external_libs/yaml-cpp)
add_subdirectory(external_libs/yaml-cpp EXCLUDE_FROM_ALL)
endif()

if (NOT TARGET yaml-cpp::yaml-cpp)
message(FATAL_ERROR "yaml-cpp was found, but no usable CMake target was provided")
endif()

# ------------------------------------------------------
# find spdlog
find_package(spdlog QUIET)
if (spdlog_FOUND AND NOT ${USE_EXTERNAL_LIBS})
message(STATUS "spdlog found")
if (NOT USE_EXTERNAL_LIBS)
find_package(spdlog QUIET)
endif()

if (spdlog_FOUND AND NOT USE_EXTERNAL_LIBS)
message(STATUS "Found system spdlog: ${spdlog_DIR}")
else()
message(STATUS "spdlog not found. use source code external_libs/spdlog...")
initialize_submodule(external_libs/spdlog)
add_subdirectory(external_libs/spdlog EXCLUDE_FROM_ALL)
message(STATUS "System spdlog not found (or bundled libraries requested); using external_libs/spdlog")
initialize_submodule(external_libs/spdlog)
add_subdirectory(external_libs/spdlog EXCLUDE_FROM_ALL)
endif()

if (TARGET spdlog::spdlog_header_only)
Expand Down
20 changes: 10 additions & 10 deletions examples/01_checkerboard_azimuthal_ani/input_params.yml
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ data:
# Love-wave source-receiver files (used when wave_type[1] is true)
src_rec_file_lv_ph: OUTPUT_FILES/src_rec_file_forward_LV_PH.csv
src_rec_file_lv_gr: OUTPUT_FILES/src_rec_file_forward_LV_GR.csv
wave_type: [True, False] # [use_rayleigh, use_love] - both required for radial anisotropy
wave_type: [True, False] # [use_rayleigh, use_love]
vel_type: [True, False] # [use_phase, use_group]
weights: [1.0, 0.5] # weights for phase and group velocity

Expand All @@ -26,9 +26,9 @@ domain:
grid_method_0:
interval: [0.02, 0.02, 0.5] # interval of lon(degree), lat(degree), z(km)
num_grid_margin: 5 # number of grid margin
grid_method_1:
lon_min_max: [119.5, 120.5]
lat_min_max: [29.5, 30.5]
grid_method_1: # only read when grid_method: 1; values below describe the same domain as grid_method_0
lon_min_max: [101.8, 104.2] # longitude range in degree
lat_min_max: [23.8, 26.2] # latitude range in degree
n_grid: [100, 100, 30] # number of grid in lon, lat, z

model:
Expand All @@ -39,13 +39,13 @@ model:

topo:
is_consider_topo: False # whether consider topography
topo_file: /path/to/topo.grd # Path to local elevation file. Only valid in topo_type: 1 and 2
topo_file: /path/to/topo.grd # Path to local elevation file. Only valid when is_consider_topo: True
wavelen_factor: 2.5

postproc:
kdensity_coe: 0 # Kdensity_coe is used to rescale the final kernel: kernel -> kernel / pow(density of kernel, Kdensity_coe).
independent_smooth_ani: False # if True, smooth azimuthal anisotropy (Gc/Gs) kernels with sigma_ani / n_inv_grid_ani
smooth_method: 0
smooth_method: 0 # 0: PDE (Gaussian) smoothing; 1: multigrid smoothing
smooth_method_0:
sigma: [0.1, 1] # smoothing lengths [horizontal (degree), vertical (km)]
sigma_ani: [0.15, 1] # same units as sigma; azimuthal anisotropy (Gc/Gs) kernels only, ignored for radial anisotropy (gamma uses sigma)
Expand All @@ -59,12 +59,12 @@ inversion:
model_para_type: 1 # 0 Isotropic, 1 Azimuthal anisotropy, 2 Radial anisotropy
use_alpha_beta_rho: false
rho_scaling: false
vpvs_ratio_range: [1.4, 2.5]
# ------- Inversion parameters --------
niter: 80 # max iterations
vpvs_ratio_range: [1.4, 2.5] # clamp range of vp/vs. Only used when use_alpha_beta_rho: True
# ------- Iteration control -------
niter: 40 # max iterations
min_derr: 0.0015 # minimum error change
# ------- optimization parameters -------
optim_method: 1 # 0: SD (adeptive step length); 1: LBFGS (line search)
optim_method: 1 # 0: SD (adaptive step length); 1: LBFGS (line search)
step_length: 0.02 # starting step length
maxshrink: 0.6 # max step length descent
# ------- Line search for LBFGS -------
Expand Down
16 changes: 8 additions & 8 deletions examples/02_checkerboard_radial_ani/input_params.yml
Original file line number Diff line number Diff line change
Expand Up @@ -26,9 +26,9 @@ domain:
grid_method_0:
interval: [0.02, 0.02, 0.2] # interval of lon(degree), lat(degree), z(km)
num_grid_margin: 5 # number of grid margin
grid_method_1:
lon_min_max: [119.5, 120.5]
lat_min_max: [29.5, 30.5]
grid_method_1: # only read when grid_method: 1; values below describe the same domain as grid_method_0
lon_min_max: [101.95, 104.05] # longitude range in degree
lat_min_max: [23.95, 26.05] # latitude range in degree
n_grid: [100, 100, 30] # number of grid in lon, lat, z

model:
Expand All @@ -39,13 +39,13 @@ model:

topo:
is_consider_topo: False # whether consider topography
topo_file: /path/to/topo.grd # Path to local elevation file. Only valid in topo_type: 1 and 2
topo_file: /path/to/topo.grd # Path to local elevation file. Only valid when is_consider_topo: True
wavelen_factor: 2.5

postproc:
kdensity_coe: 0 # Kdensity_coe is used to rescale the final kernel: kernel -> kernel / pow(density of kernel, Kdensity_coe).
independent_smooth_ani: False # if True, smooth azimuthal anisotropy (Gc/Gs) kernels with sigma_ani / n_inv_grid_ani
smooth_method: 0
smooth_method: 0 # 0: PDE (Gaussian) smoothing; 1: multigrid smoothing
smooth_method_0:
sigma: [0.1, 2] # smoothing lengths [horizontal (degree), vertical (km)]
sigma_ani: [0.1, 2] # same units as sigma; azimuthal anisotropy (Gc/Gs) kernels only, ignored for radial anisotropy (gamma uses sigma)
Expand All @@ -59,12 +59,12 @@ inversion:
model_para_type: 2 # 0 Isotropic, 1 Azimuthal anisotropy, 2 Radial anisotropy
use_alpha_beta_rho: false
rho_scaling: false
vpvs_ratio_range: [1.4, 2.5]
# ------- Inversion parameters --------
vpvs_ratio_range: [1.4, 2.5] # clamp range of vp/vs. Only used when use_alpha_beta_rho: True
# ------- Iteration control -------
niter: 40 # max iterations
min_derr: 0.0015 # minimum error change
# ------- optimization parameters -------
optim_method: 1 # 0: SD (adeptive step length); 1: LBFGS (line search)
optim_method: 1 # 0: SD (adaptive step length); 1: LBFGS (line search)
step_length: 0.02 # starting step length
maxshrink: 0.6 # max step length descent
# ------- Line search for LBFGS -------
Expand Down
21 changes: 13 additions & 8 deletions examples/04_hawaii_topo/input_params.yml
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,10 @@ domain:
grid_method_0:
interval: [0.01, 0.01, 0.5] # interval of lon(degree), lat(degree), z(km)
num_grid_margin: 5 # number of grid margin
grid_method_1: # only read when grid_method: 1; values below describe the same domain as grid_method_0
lon_min_max: [-0.52, 0.66] # longitude range in degree (rotated frame, see rotatation.ipynb)
lat_min_max: [-0.47, 0.48] # latitude range in degree (rotated frame)
n_grid: [120, 100, 31] # number of grid in lon, lat, z

model:
# ------- Initial model --------
Expand All @@ -31,16 +35,16 @@ model:

topo:
is_consider_topo: True # whether consider topography
topo_file: hawaii_rotated.nc # Path to local elevation file. Only valid in topo_type: 1 and 2
topo_file: hawaii_rotated.nc # Path to local elevation file. Only valid when is_consider_topo: True
wavelen_factor: 2.5

postproc:
kdensity_coe: 1 # Kdensity_coe is used to rescale the final kernel: kernel -> kernel / pow(density of kernel, Kdensity_coe).
kdensity_coe: 0.0 # Kdensity_coe is used to rescale the final kernel: kernel -> kernel / pow(density of kernel, Kdensity_coe).
independent_smooth_ani: False # if True, smooth azimuthal anisotropy (Gc/Gs) kernels with sigma_ani / n_inv_grid_ani
smooth_method: 1 # 0: PDE (Gaussian) smoothing; 1: multigrid smoothing
smooth_method_0:
sigma: [0.1, 2] # smoothing lengths [horizontal (degree), vertical (km)]
sigma_ani: [0.1, 2] # same units as sigma; azimuthal anisotropy (Gc/Gs) kernels only, ignored for radial anisotropy (gamma uses sigma)
sigma: [0.12, 1.5] # smoothing lengths [horizontal (degree), vertical (km)]
sigma_ani: [0.12, 1.5] # same units as sigma; azimuthal anisotropy (Gc/Gs) kernels only, ignored for radial anisotropy (gamma uses sigma)
smooth_method_1:
n_inv_components: 5
n_inv_grid: [8, 9, 10]
Expand All @@ -49,13 +53,14 @@ postproc:
inversion:
# ------- Inversion parameters --------
model_para_type: 0 # 0 Isotropic, 1 Azimuthal anisotropy, 2 Radial anisotropy
use_alpha_beta_rho: true
rho_scaling: true
# ------- Inversion parameters --------
use_alpha_beta_rho: false
rho_scaling: false
vpvs_ratio_range: [1.4, 2.5] # clamp range of vp/vs. Only used when use_alpha_beta_rho: True
# ------- Iteration control -------
niter: 40 # max iterations
min_derr: 0.0001 # minimum error change
# ------- optimization parameters -------
optim_method: 1 # 0: SD (adeptive step length); 1: LBFGS (line search)
optim_method: 0 # 0: SD (adaptive step length); 1: LBFGS (line search)
step_length: 0.02 # starting step length
maxshrink: 0.6 # max step length descent
# ------- Line search for LBFGS -------
Expand Down
55 changes: 9 additions & 46 deletions examples/04_hawaii_topo/plot_model.ipynb

Large diffs are not rendered by default.

55 changes: 11 additions & 44 deletions include/argparser.h
Original file line number Diff line number Diff line change
Expand Up @@ -307,19 +307,14 @@ inline RotateTopoArgs argparse_rotate_topo(int argc, char* argv[]) {
};
}

// Model fields that surfatt_rotate_model is allowed to export. Every entry is a
// scalar field, i.e. invariant under a rotation of the horizontal frame, so only
// the grid coordinates have to be transformed and the values travel unchanged.
inline const std::array<const char*, 7> rotate_model_fields = {
"vs", "vsv", "vsh", "vp", "rho", "zeta", "gamma"
// Model fields written by ModelGrid::write(), plus gamma used by model_iter.h5.
// SURFATT_rotate_model passes scalar fields through unchanged and rotates the
// frame-dependent gc/gs/theta fields together with the grid coordinates.
inline const std::array<const char*, 11> rotate_model_fields = {
"vs", "vp", "rho", "gc", "gs", "g0", "theta",
"vsv", "vsh", "zeta", "gamma"
};

// Fields that are NOT invariant under a frame rotation and must never be exported
// as a plain pass-through: gc/gs are the cos(2*theta)/sin(2*theta) components of
// the azimuthal anisotropy tensor, and theta is the fast-axis azimuth itself (it
// is rotated, but only on the automatic g0/theta columns).
inline const std::array<const char*, 3> rotate_model_blocked_fields = {"gc", "gs", "theta"};

// Comma-separated list of rotate_model_fields, for help text and error messages.
inline std::string rotate_model_fields_str() {
std::string s;
Expand Down Expand Up @@ -359,15 +354,15 @@ inline RotateModelArgs argparse_rotate_model(int argc, char* argv[]) {
" -c clat/clon Centre of rotation (lat/lon)\n\n"
"optional arguments:\n"
" -a angle Rotation angle in degrees (default: 0)\n"
" -k keyname Export only this scalar dataset as the model column\n"
" -k keyname Export only this dataset as the model column\n"
" One of: " << rotate_model_fields_str() << "\n"
" optionally with a model_/grad_ prefix and an _NNN\n"
" iteration suffix, e.g. model_vs_042 in model_iter.h5\n"
" The csv column is named after the bare field, so\n"
" model_vs_042 is written as a column named vs\n"
" Without -k, vs is exported together with whichever\n"
" anisotropy fields the file carries: g0/theta for\n"
" azimuthal, vsv/zeta for radial\n"
" Without -k, every recognised model field present\n"
" in the file is exported. gc/gs/theta are rotated\n"
" together with the coordinate frame.\n"
" -h Print help message\n";
std::exit(0);
}
Expand All @@ -376,34 +371,7 @@ inline RotateModelArgs argparse_rotate_model(int argc, char* argv[]) {
out.outfname = al.require("-o");
if (auto v = al.get("-a")) out.angle = std::stod(*v);
if (auto v = al.get("-c")) out.center = parse_2double(*v);
if (auto v = al.get("-k")) {
const std::string base = rotate_model_base_field(*v);
const auto matches = [&base](const char* k) { return base == k; };
// Reported directly instead of thrown: an uncaught exception under mpirun
// buries the message in an abort trace, and this is a routine user error.
if (std::any_of(rotate_model_blocked_fields.begin(),
rotate_model_blocked_fields.end(), matches)) {
std::cerr << "Error: -k " << *v << " is not supported: \"" << base
<< "\" is not invariant under a rotation of the frame.\n";
if (base == "theta")
std::cerr << "theta is the fast-axis azimuth; it is rotated together with the grid "
"and is\nexported next to g0 when -k is omitted and both are in the "
"file.\n";
else
std::cerr << "gc/gs are the cos(2*theta)/sin(2*theta) components of the azimuthal "
"anisotropy\ntensor; the derived g0/theta pair is exported when -k is "
"omitted.\n";
std::exit(1);
}
if (!std::any_of(rotate_model_fields.begin(), rotate_model_fields.end(), matches)) {
std::cerr << "Error: unknown dataset \"" << *v << "\" for -k.\n"
<< "Allowed fields: " << rotate_model_fields_str() << "\n"
<< "optionally with a model_/grad_ prefix and an _NNN iteration suffix, "
"e.g. model_vs_042\n";
std::exit(1);
}
out.key = *v;
}
if (auto v = al.get("-k")) out.key = *v;
return out;
}

Expand Down Expand Up @@ -432,4 +400,3 @@ inline Tomo2DArgs argparse_tomo2d(int argc, char* argv[]) {
if (auto v = al.get("-m")) out.hmarg = std::stod(*v);
return out;
}

Loading
Loading