From bc703906ce9451462adcc34dffadb089d9747e7e Mon Sep 17 00:00:00 2001 From: xumi1993 Date: Tue, 18 Aug 2026 21:46:33 +0800 Subject: [PATCH 1/6] Update project version to 2.1.2 and add kernel density storage functionality --- CMakeLists.txt | 2 +- include/inversion.h | 1 + src/inversion.cpp | 28 ++++++++++++++++++++++++++++ src/xdmf.cpp | 3 +++ 4 files changed, 33 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index f092c0d..5ae837e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ cmake_minimum_required(VERSION 3.18) -project( SurfATT VERSION 2.1.1 LANGUAGES CXX Fortran) +project( SurfATT VERSION 2.1.2 LANGUAGES CXX Fortran) include(GNUInstallDirs) diff --git a/include/inversion.h b/include/inversion.h index 235af41..cc03f69 100644 --- a/include/inversion.h +++ b/include/inversion.h @@ -59,6 +59,7 @@ class Inversion { void grad_normalization(FieldVec &grads); void store_model(); void store_gradient(); + void store_kernel_density(); bool check_convergence(); void model_update(FieldVec &dir); bool line_search(); diff --git a/src/inversion.cpp b/src/inversion.cpp index 22b0721..0279c5f 100644 --- a/src/inversion.cpp +++ b/src/inversion.cpp @@ -199,6 +199,10 @@ void Inversion::run_inversion() { write_obj_line(); + if (IP.postproc().is_kden) { + store_kernel_density(); + } + if (IP.output().output_in_process_model || IP.inversion().optim_method == OPTIM_LBFGS) { store_gradient(); } @@ -496,6 +500,30 @@ void Inversion::store_gradient() { xdmf::write_model_iter(xdmf_fname_, iter_, iter_); } +// Save the total kernel density used during the current inversion iteration. +// Each active data type contributes according to its phase/group data weight. +// Dataset names follow the model/gradient convention: kernel_density_{N}. +void Inversion::store_kernel_density() { + auto &IP = InputParams::IP(); + auto &dcp = Decomposer::DCP(); + auto &mpi = Parallel::mpi(); + + Tensor3r density_loc(dcp.loc_nx(), dcp.loc_ny(), ngrid_k); + density_loc.setZero(); + for (auto [wt, tp] : IP.data().active_data) { + const int itype = static_cast(tp); + density_loc += SurfGrid::SG(wt, tp).ker_den_loc * IP.data().weights[itype]; + } + + // collect_data is collective: all ranks participate, then only the main + // rank owns and writes the assembled global volume. + Tensor3r density_all = dcp.collect_data(density_loc.data()); + if (!mpi.is_main()) return; + + H5IO f(db_fname, H5IO::RDWR); + f.write_tensor(fmt::format("kernel_density_{:03d}", iter_), density_all); +} + void Inversion::steepest_descent() { auto &IP = InputParams::IP(); auto &logger = ATTLogger::logger(); diff --git a/src/xdmf.cpp b/src/xdmf.cpp index 28065ae..41656be 100644 --- a/src/xdmf.cpp +++ b/src/xdmf.cpp @@ -93,6 +93,9 @@ void write_model_iter(const std::string &xdmf_path, int iter, out << attr("grad_gc", "grad_gc" + sfx) << attr("grad_gs", "grad_gs" + sfx); } + if (IP.postproc().is_kden) { + out << attr("kernel_density", "kernel_density" + sfx); + } } out << " \n"; } From 6c6b5d3da0bee4f7926c02b6086aa6a2bf33aab2 Mon Sep 17 00:00:00 2001 From: xumi1993 Date: Fri, 11 Sep 2026 17:14:35 +0800 Subject: [PATCH 2/6] Enhance radial anisotropy handling: update Vs calculation, add Vsh support, and refine kernel computations --- include/utils.h | 17 +++++++++- src/inversion.cpp | 34 ++++++++------------ src/model_grid.cpp | 73 +++++++++++++++++++++++++------------------ src/postproc.cpp | 2 +- src/preproc.cpp | 31 ++++++++++++++++-- src/surf_grid.cpp | 18 +++++++---- src/surfker/rlker.cpp | 13 +++++--- 7 files changed, 121 insertions(+), 67 deletions(-) diff --git a/include/utils.h b/include/utils.h index 4737ec4..3a7a3f1 100644 --- a/include/utils.h +++ b/include/utils.h @@ -680,7 +680,7 @@ inline std::vector fmt_col(const real_t* data, int n, int prec = 6) // Converts between (vsv, vsh) and (Vs, zeta) parameterizations. // // Definitions: -// Vs = sqrt((2*vsv + vsh) / 3) — RMS S-wave velocity +// Vs = sqrt((2*vsv^2 + vsh^2) / 3) — RMS S-wave velocity // zeta = vsh^2 / vsv^2 — anisotropy ratio // // Supports scalar, vector (Eigen::VectorX), and tensor (Eigen::Tensor<3>) inputs. @@ -691,6 +691,21 @@ inline real_t vsvvsh2vs(real_t vsv, real_t vsh) { return std::sqrt((2 * vsv* vsv + vsh * vsh) / 3); } +// Absolute partial derivatives at independent (Vsv, Vsh), including the +// shared empirical Vp(Vs) and rho(Vp(Vs)). The dispersion solver's beta is +// Vsv for Rayleigh and Vsh for Love. Convert to log coordinates with +// K_lnVsv = Vsv*K_vsv + Vsh*K_vsh, K_lngamma = Vsh*K_vsh. +inline std::pair radial_material_kernels( + real_t vsv, real_t vsh, real_t k_beta, real_t k_vp, + real_t k_rho, bool love +) { + const real_t vs = vsvvsh2vs(vsv, vsh); + const real_t shared = dalpha_dbeta(vs) * + (k_vp + k_rho * drho_dalpha(vs2vp(vs))); + return {(love ? _0_CR : k_beta) + shared * 2 * vsv / (3 * vs), + (love ? k_beta : _0_CR) + shared * vsh / (3 * vs)}; +} + // Scalar version: compute zeta from vsv and vsh inline real_t vsvvsh2zeta(real_t vsv, real_t vsh) { return (vsh * vsh) / (vsv * vsv); diff --git a/src/inversion.cpp b/src/inversion.cpp index 0279c5f..44427c9 100644 --- a/src/inversion.cpp +++ b/src/inversion.cpp @@ -329,14 +329,9 @@ void Inversion::accumulate_smoothed_gradient( }; if (model_para_type == MODEL_RADIAL_ANI) { - if (wt == WaveType::RL) { - for (int ipara = 0; ipara < 3; ++ipara) { - accumulate_grad(ipara, true); - } - } else if (wt == WaveType::LV) { - // Keep previous behavior: gamma term is always accumulated for Love in radial anisotropy. - accumulate_grad(NPARAMS - 1, false); - } + // Both wave types contribute independent Vsv and Vsh partials. + accumulate_grad(0, true); + accumulate_grad(5, false); // also needed by the existing radial conversion return; } @@ -380,15 +375,9 @@ real_t Inversion::run_forward_adjoint(const bool is_calc_adj) { if (IP.inversion().optim_method == OPTIM_LBFGS) { const real_t wt_val = IP.data().weights[itype]; if (IP.inversion().model_para_type == MODEL_RADIAL_ANI) { - // For radial anisotropy, Love kernel lives in ker_loc[0] and contributes - // to the gamma (index 5) direction; Rayleigh only updates vs/vp/rho. - if (wt == WaveType::RL) { - for (int ipara = 0; ipara < 3; ++ipara) - if (is_active_param[ipara]) - ker_curr_[ipara] = ker_curr_[ipara] + sg.ker_loc[ipara] * wt_val; - } else { - ker_curr_[5] = ker_curr_[5] + sg.ker_loc[0] * wt_val; - } + // Physical partials in independent (Vsv, Vsh), for both waves. + ker_curr_[0] += sg.ker_loc[0] * wt_val; + ker_curr_[5] += sg.ker_loc[5] * wt_val; } else { for (int ipara = 0; ipara < NPARAMS; ++ipara) if (is_active_param[ipara]) @@ -653,10 +642,6 @@ void Inversion::model_update(FieldVec &dir) { mg.vp3d_loc = mg.vp3d_loc * (1 - alpha_ * dir[1]); mg.rho3d_loc = mg.rho3d_loc * (1 - alpha_ * dir[2]); alpha_clamp(); - } else { - // Empirical scaling: vs → vp → rho via Brocher (2005) - mg.vp3d_loc = vs2vp(mg.vs3d_loc); - mg.rho3d_loc = vp2rho(mg.vp3d_loc); } if (IP.inversion().model_para_type == MODEL_AZI_ANI) { mg.gc3d_loc = mg.gc3d_loc - alpha_ * dir[3]; @@ -666,6 +651,13 @@ void Inversion::model_update(FieldVec &dir) { mg.gamma3d_loc = mg.gamma3d_loc * (1 - alpha_ * dir[5]); mg.vsh3d_loc = mg.vs3d_loc * mg.gamma3d_loc; } + if (!IP.inversion().use_alpha_beta_rho) { + // Recompute only after both Vsv and Vsh have been updated. + const Tensor3r mean_vs = IP.inversion().model_para_type == MODEL_RADIAL_ANI + ? vsvvsh2vs(mg.vs3d_loc, mg.vsh3d_loc) : mg.vs3d_loc; + mg.vp3d_loc = vs2vp(mean_vs); + mg.rho3d_loc = vp2rho(mg.vp3d_loc); + } mpi.barrier(); } diff --git a/src/model_grid.cpp b/src/model_grid.cpp index c99949e..888fd3b 100644 --- a/src/model_grid.cpp +++ b/src/model_grid.cpp @@ -412,34 +412,6 @@ void ModelGrid::load_3d_model() { ); } - // --- vp (optional, fallback: empirical vs2vp) --- - if (f.exists("vp")) { - interpolate_or_copy("vp", vp3d); - logger.Info( - same_grid ? "Loaded 'vp' from HDF5 file." : - "Loaded and interpolated 'vp' onto current model grid.", - MODULE_GRID - ); - } else { - logger.Info("'vp' not found in HDF5 file, computing from empirical vs2vp.", MODULE_GRID); - for (hsize_t i = 0; i < expect_n; ++i) - vp3d[i] = vs2vp(vs3d[i]); - } - - // --- rho (optional, fallback: empirical vp2rho) --- - if (f.exists("rho")) { - interpolate_or_copy("rho", rho3d); - logger.Info( - same_grid ? "Loaded 'rho' from HDF5 file." : - "Loaded and interpolated 'rho' onto current model grid.", - MODULE_GRID - ); - } else { - logger.Info("'rho' not found in HDF5 file, computing from empirical vp2rho.", MODULE_GRID); - for (hsize_t i = 0; i < expect_n; ++i) - rho3d[i] = vp2rho(vp3d[i]); - } - // --- gc / gs (optional, only if model_para_type is MODEL_AZI_ANI) --- if (IP.inversion().model_para_type == MODEL_AZI_ANI) { if (f.exists("gc")) { @@ -490,6 +462,45 @@ void ModelGrid::load_3d_model() { } } } + + // Initialize Vp and density once, after both shear components are ready. + // Radial inversion derives these fields from mean Vs; other modes retain + // the optional input fields and their empirical fallbacks. + if (IP.inversion().model_para_type == MODEL_RADIAL_ANI) { + logger.Info("Computing shared Vp and density from mean Vs for the radial model.", MODULE_GRID); + for (hsize_t i = 0; i < expect_n; ++i) { + vp3d[i] = vs2vp(vsvvsh2vs(vs3d[i], vsh3d[i])); + rho3d[i] = vp2rho(vp3d[i]); + } + } else { + // --- vp (optional, fallback: empirical vs2vp) --- + if (f.exists("vp")) { + interpolate_or_copy("vp", vp3d); + logger.Info( + same_grid ? "Loaded 'vp' from HDF5 file." : + "Loaded and interpolated 'vp' onto current model grid.", + MODULE_GRID + ); + } else { + logger.Info("'vp' not found in HDF5 file, computing from empirical vs2vp.", MODULE_GRID); + for (hsize_t i = 0; i < expect_n; ++i) + vp3d[i] = vs2vp(vs3d[i]); + } + + // --- rho (optional, fallback: empirical vp2rho) --- + if (f.exists("rho")) { + interpolate_or_copy("rho", rho3d); + logger.Info( + same_grid ? "Loaded 'rho' from HDF5 file." : + "Loaded and interpolated 'rho' onto current model grid.", + MODULE_GRID + ); + } else { + logger.Info("'rho' not found in HDF5 file, computing from empirical vp2rho.", MODULE_GRID); + for (hsize_t i = 0; i < expect_n; ++i) + rho3d[i] = vp2rho(vp3d[i]); + } + } } catch (const std::exception &e) { logger.Error(fmt::format( "ModelGrid: failed to load 3D model from HDF5 file '{}': {}", @@ -738,10 +749,10 @@ void ModelGrid::add_radial_aniso_perturbation( // vsh = vsv * sqrt(zeta) auto [vsv_new, vsh_new] = recover_anisotropy(Vs_new, zeta_new); - // Update vs3d to store Vs (for consistency with load_3d_model) + // vs3d stores Vsv; Vp and density use the mean Vs. vs3d[idx] = vsv_new; vsh3d[idx] = vsh_new; - vp3d[idx] = vs2vp(vsv_new); + vp3d[idx] = vs2vp(Vs_new); rho3d[idx] = vp2rho(vp3d[idx]); } } @@ -816,7 +827,7 @@ void ModelGrid::write(const std::string &subname) { f.write_volume("vsv", vs3d, ngrid_i, ngrid_j, ngrid_k); f.write_volume("vsh", vsh3d, ngrid_i, ngrid_j, ngrid_k); - // Compute and write Vs = sqrt((2*vsv + vsh)/3) and zeta = vsh^2/vsv^2 + // Compute and write Vs = sqrt((2*vsv^2 + vsh^2)/3) and zeta = vsh^2/vsv^2 std::vector Vs(nelem); std::vector zeta(nelem); for (int i = 0; i < nelem; ++i) { diff --git a/src/postproc.cpp b/src/postproc.cpp index b25d8a3..147b2ea 100644 --- a/src/postproc.cpp +++ b/src/postproc.cpp @@ -518,7 +518,7 @@ void postproc::kernel_precondition(SurfGrid& sg) { ); // Precondition the kernels by multiplying with the reference model parameters at each surface grid point - int nker = NPARAMS - 1; + int nker = NPARAMS; for (int iparam = 0; iparam < nker; ++iparam) { if (sg.is_active_ker(iparam)) { sg.ker_loc[iparam] = sg.ker_loc[iparam] * hess_inv; diff --git a/src/preproc.cpp b/src/preproc.cpp index 5a3fdad..929696c 100644 --- a/src/preproc.cpp +++ b/src/preproc.cpp @@ -212,6 +212,10 @@ void preproc::combine_kernels(SurfGrid& sg) { // vs kernel — always allocated (index 0) sg.ker_loc[0] = Tensor3r(dcp.loc_nx(), dcp.loc_ny(), ngrid_k); sg.ker_loc[0].setZero(); + if (IP.inversion().model_para_type == MODEL_RADIAL_ANI) { + sg.ker_loc[5] = Tensor3r(dcp.loc_nx(), dcp.loc_ny(), ngrid_k); + sg.ker_loc[5].setZero(); + } if (IP.inversion().use_alpha_beta_rho && sg.wave_type() == WaveType::RL) { // vp (1) and rho (2) — only when parametrised independently sg.ker_loc[1] = Tensor3r(dcp.loc_nx(), dcp.loc_ny(), ngrid_k); @@ -251,7 +255,31 @@ void preproc::combine_kernels(SurfGrid& sg) { } // Isotropic parameter kernels — gated by use_alpha_beta_rho only - if (IP.inversion().use_alpha_beta_rho && sg.wave_type() == WaveType::RL) { + if (IP.inversion().model_para_type == MODEL_RADIAL_ANI) { + for (int ix = 0; ix < dcp.loc_nx(); ++ix) { + for (int iy = 0; iy < dcp.loc_ny(); ++iy) { + const int gx = dcp.loc_I_start() + ix; + const int gy = dcp.loc_J_start() + iy; + for (int k = 0; k < ngrid_k; ++k) { + const real_t vsv = mg.vs3d_loc(ix, iy, k); + const real_t vsh = mg.vsh3d_loc(ix, iy, k); + const auto [ka, kb] = radial_material_kernels( + vsv, vsh, sg.sen_vs_loc(ix, iy, k, iper), + sg.sen_vp_loc(ix, iy, k, iper), + sg.sen_rho_loc(ix, iy, k, iper), + sg.wave_type() == WaveType::LV); + sg.ker_loc[0](ix, iy, k) -= adj_tt(gx, gy) * ka; + sg.ker_loc[5](ix, iy, k) -= adj_tt(gx, gy) * kb; + if (IP.postproc().is_kden) { + // Preserve the wave-specific shear-speed coverage scale, + // including the required empirical-material chain rule. + sg.ker_den_loc(ix, iy, k) -= adj_den(gx, gy) * + (sg.wave_type() == WaveType::LV ? kb : ka); + } + } + } + } + } else if (IP.inversion().use_alpha_beta_rho && sg.wave_type() == WaveType::RL) { for (int ix = 0; ix < dcp.loc_nx(); ++ix) { for (int iy = 0; iy < dcp.loc_ny(); ++iy) { const int iglob_x = dcp.loc_I_start() + ix; @@ -358,4 +386,3 @@ void preproc::combine_kernels(SurfGrid& sg) { mpi.barrier(); } - diff --git a/src/surf_grid.cpp b/src/surf_grid.cpp index 1267bd0..eb3e918 100644 --- a/src/surf_grid.cpp +++ b/src/surf_grid.cpp @@ -82,10 +82,12 @@ SurfGrid::SurfGrid(WaveType wt, SurfType vt){ void SurfGrid::setup_active_kernels() { auto& IP = InputParams::IP(); - active_kernels_ = std::vector(NPARAMS-1, false); - active_kernels_[0] = true; // vp + active_kernels_ = std::vector(NPARAMS, false); + active_kernels_[0] = true; // vs (Vsv for radial models) + if (IP.inversion().model_para_type == MODEL_RADIAL_ANI) + active_kernels_[5] = true; // independent Vsh partial, for both waves if (IP.inversion().use_alpha_beta_rho && wt_ == WaveType::RL) { - active_kernels_[1] = true; // vs + active_kernels_[1] = true; // vp active_kernels_[2] = true; // rho } if (IP.inversion().model_para_type == MODEL_AZI_ANI && wt_ == WaveType::RL) { @@ -298,7 +300,10 @@ void SurfGrid::compute_dispersion_kernel() { vp1d(k) = mg.vp3d_loc(ix, iy, k); rho1d(k) = mg.rho3d_loc(ix, iy, k); } else { - vp1d(k) = vs2vp(vs1d(k)); + const real_t mean_vs = IP.inversion().model_para_type == MODEL_RADIAL_ANI + ? vsvvsh2vs(mg.vs3d_loc(ix, iy, k), mg.vsh3d_loc(ix, iy, k)) + : vs1d(k); + vp1d(k) = vs2vp(mean_vs); rho1d(k) = vp2rho(vp1d(k)); } } @@ -317,10 +322,11 @@ void SurfGrid::compute_dispersion_kernel() { vs_block = kernels.sen_vs.transpose(); if (wt_ == WaveType::RL){ Eigen::Map vp_block(sen_vp_loc.data() + id0, ngrid_k, nperiod_); - Eigen::Map rho_block(sen_rho_loc.data() + id0, ngrid_k, nperiod_); vp_block = kernels.sen_vp.transpose(); - rho_block = kernels.sen_rho.transpose(); } + // Love also has a density kernel; radial empirical scaling needs it. + Eigen::Map rho_block(sen_rho_loc.data() + id0, ngrid_k, nperiod_); + rho_block = kernels.sen_rho.transpose(); if (IP.inversion().model_para_type == MODEL_AZI_ANI) { Eigen::Map gc_block(sen_gc_loc.data() + id0, ngrid_k, nperiod_); Eigen::Map gs_block(sen_gs_loc.data() + id0, ngrid_k, nperiod_); diff --git a/src/surfker/rlker.cpp b/src/surfker/rlker.cpp index 02381f5..acb6d8a 100644 --- a/src/surfker/rlker.cpp +++ b/src/surfker/rlker.cpp @@ -790,7 +790,7 @@ static void energy_func(double om, double wvno, RayleighWork& s) double facav = rho * av * av * DUZDUZ / wvno2; s.dcda[m] = facah + facav; double facr = -0.5 * c * c * (URUR + UZUZ); - s.dcdr[m] = 0.5 * (av * facav + ah * facah) + facr * rho; + s.dcdr[m] = 0.5 * (facav + facah) + facr * rho; s.dcdb[m] = 0.0; s.dcdgc[m] = 0.0; s.dcdgs[m] = 0.0; @@ -819,7 +819,7 @@ static void energy_func(double om, double wvno, RayleighWork& s) s.dcda[m] = facah + facav; s.dcdb[m] = facbv + facbh; double facr = -0.5 * c * c * (URUR + UZUZ); - s.dcdr[m] = 0.5 * (av * facav + ah * facah + bv * facbv) + facr * rho; + s.dcdr[m] = 0.5 * (facav + facah + facbv) + facr * rho; s.dcdgc[m] = rho*bv*bv*(UZUZ + 2.0*UZDUR/wvno + DURDUR/wvno2) + rho*bv*bv*URUR; s.dcdgs[m] = s.dcdgc[m]; } @@ -833,9 +833,12 @@ static void energy_func(double om, double wvno, RayleighWork& s) double inv_ug_s0 = 1.0 / (s.ugr * s.sumi0); for (int m = 0; m < mmax; m++) { - s.dcda[m] *= inv_ug_s0; - s.dcdb[m] *= inv_ug_s0; - s.dcdr[m] *= inv_ug_s0; + // facah/facav/facbv contain rho * velocity^2: dividing only by + // U*I0 gives logarithmic velocity derivatives. The public kernels + // are absolute derivatives dc/dVp, dc/dVs and dc/drho. + s.dcda[m] *= inv_ug_s0 / s.za[m]; + s.dcdb[m] *= (s.iwat[m] == 1) ? 0.0 : inv_ug_s0 / s.zb[m]; + s.dcdr[m] *= inv_ug_s0 / s.zrho[m]; s.dcdgc[m] *= inv_ug_s0; s.dcdgs[m] *= inv_ug_s0; } From 41797967ad3ecfeb5ba026018c6cb254311a7fe7 Mon Sep 17 00:00:00 2001 From: xumi1993 Date: Fri, 11 Sep 2026 20:41:07 +0800 Subject: [PATCH 3/6] Revert "Enhance radial anisotropy handling: update Vs calculation, add Vsh support, and refine kernel computations" This reverts commit 6c6b5d3da0bee4f7926c02b6086aa6a2bf33aab2. --- include/utils.h | 17 +--------- src/inversion.cpp | 34 ++++++++++++-------- src/model_grid.cpp | 73 ++++++++++++++++++------------------------- src/postproc.cpp | 2 +- src/preproc.cpp | 31 ++---------------- src/surf_grid.cpp | 18 ++++------- src/surfker/rlker.cpp | 13 +++----- 7 files changed, 67 insertions(+), 121 deletions(-) diff --git a/include/utils.h b/include/utils.h index 3a7a3f1..4737ec4 100644 --- a/include/utils.h +++ b/include/utils.h @@ -680,7 +680,7 @@ inline std::vector fmt_col(const real_t* data, int n, int prec = 6) // Converts between (vsv, vsh) and (Vs, zeta) parameterizations. // // Definitions: -// Vs = sqrt((2*vsv^2 + vsh^2) / 3) — RMS S-wave velocity +// Vs = sqrt((2*vsv + vsh) / 3) — RMS S-wave velocity // zeta = vsh^2 / vsv^2 — anisotropy ratio // // Supports scalar, vector (Eigen::VectorX), and tensor (Eigen::Tensor<3>) inputs. @@ -691,21 +691,6 @@ inline real_t vsvvsh2vs(real_t vsv, real_t vsh) { return std::sqrt((2 * vsv* vsv + vsh * vsh) / 3); } -// Absolute partial derivatives at independent (Vsv, Vsh), including the -// shared empirical Vp(Vs) and rho(Vp(Vs)). The dispersion solver's beta is -// Vsv for Rayleigh and Vsh for Love. Convert to log coordinates with -// K_lnVsv = Vsv*K_vsv + Vsh*K_vsh, K_lngamma = Vsh*K_vsh. -inline std::pair radial_material_kernels( - real_t vsv, real_t vsh, real_t k_beta, real_t k_vp, - real_t k_rho, bool love -) { - const real_t vs = vsvvsh2vs(vsv, vsh); - const real_t shared = dalpha_dbeta(vs) * - (k_vp + k_rho * drho_dalpha(vs2vp(vs))); - return {(love ? _0_CR : k_beta) + shared * 2 * vsv / (3 * vs), - (love ? k_beta : _0_CR) + shared * vsh / (3 * vs)}; -} - // Scalar version: compute zeta from vsv and vsh inline real_t vsvvsh2zeta(real_t vsv, real_t vsh) { return (vsh * vsh) / (vsv * vsv); diff --git a/src/inversion.cpp b/src/inversion.cpp index 44427c9..0279c5f 100644 --- a/src/inversion.cpp +++ b/src/inversion.cpp @@ -329,9 +329,14 @@ void Inversion::accumulate_smoothed_gradient( }; if (model_para_type == MODEL_RADIAL_ANI) { - // Both wave types contribute independent Vsv and Vsh partials. - accumulate_grad(0, true); - accumulate_grad(5, false); // also needed by the existing radial conversion + if (wt == WaveType::RL) { + for (int ipara = 0; ipara < 3; ++ipara) { + accumulate_grad(ipara, true); + } + } else if (wt == WaveType::LV) { + // Keep previous behavior: gamma term is always accumulated for Love in radial anisotropy. + accumulate_grad(NPARAMS - 1, false); + } return; } @@ -375,9 +380,15 @@ real_t Inversion::run_forward_adjoint(const bool is_calc_adj) { if (IP.inversion().optim_method == OPTIM_LBFGS) { const real_t wt_val = IP.data().weights[itype]; if (IP.inversion().model_para_type == MODEL_RADIAL_ANI) { - // Physical partials in independent (Vsv, Vsh), for both waves. - ker_curr_[0] += sg.ker_loc[0] * wt_val; - ker_curr_[5] += sg.ker_loc[5] * wt_val; + // For radial anisotropy, Love kernel lives in ker_loc[0] and contributes + // to the gamma (index 5) direction; Rayleigh only updates vs/vp/rho. + if (wt == WaveType::RL) { + for (int ipara = 0; ipara < 3; ++ipara) + if (is_active_param[ipara]) + ker_curr_[ipara] = ker_curr_[ipara] + sg.ker_loc[ipara] * wt_val; + } else { + ker_curr_[5] = ker_curr_[5] + sg.ker_loc[0] * wt_val; + } } else { for (int ipara = 0; ipara < NPARAMS; ++ipara) if (is_active_param[ipara]) @@ -642,6 +653,10 @@ void Inversion::model_update(FieldVec &dir) { mg.vp3d_loc = mg.vp3d_loc * (1 - alpha_ * dir[1]); mg.rho3d_loc = mg.rho3d_loc * (1 - alpha_ * dir[2]); alpha_clamp(); + } else { + // Empirical scaling: vs → vp → rho via Brocher (2005) + mg.vp3d_loc = vs2vp(mg.vs3d_loc); + mg.rho3d_loc = vp2rho(mg.vp3d_loc); } if (IP.inversion().model_para_type == MODEL_AZI_ANI) { mg.gc3d_loc = mg.gc3d_loc - alpha_ * dir[3]; @@ -651,13 +666,6 @@ void Inversion::model_update(FieldVec &dir) { mg.gamma3d_loc = mg.gamma3d_loc * (1 - alpha_ * dir[5]); mg.vsh3d_loc = mg.vs3d_loc * mg.gamma3d_loc; } - if (!IP.inversion().use_alpha_beta_rho) { - // Recompute only after both Vsv and Vsh have been updated. - const Tensor3r mean_vs = IP.inversion().model_para_type == MODEL_RADIAL_ANI - ? vsvvsh2vs(mg.vs3d_loc, mg.vsh3d_loc) : mg.vs3d_loc; - mg.vp3d_loc = vs2vp(mean_vs); - mg.rho3d_loc = vp2rho(mg.vp3d_loc); - } mpi.barrier(); } diff --git a/src/model_grid.cpp b/src/model_grid.cpp index 888fd3b..c99949e 100644 --- a/src/model_grid.cpp +++ b/src/model_grid.cpp @@ -412,6 +412,34 @@ void ModelGrid::load_3d_model() { ); } + // --- vp (optional, fallback: empirical vs2vp) --- + if (f.exists("vp")) { + interpolate_or_copy("vp", vp3d); + logger.Info( + same_grid ? "Loaded 'vp' from HDF5 file." : + "Loaded and interpolated 'vp' onto current model grid.", + MODULE_GRID + ); + } else { + logger.Info("'vp' not found in HDF5 file, computing from empirical vs2vp.", MODULE_GRID); + for (hsize_t i = 0; i < expect_n; ++i) + vp3d[i] = vs2vp(vs3d[i]); + } + + // --- rho (optional, fallback: empirical vp2rho) --- + if (f.exists("rho")) { + interpolate_or_copy("rho", rho3d); + logger.Info( + same_grid ? "Loaded 'rho' from HDF5 file." : + "Loaded and interpolated 'rho' onto current model grid.", + MODULE_GRID + ); + } else { + logger.Info("'rho' not found in HDF5 file, computing from empirical vp2rho.", MODULE_GRID); + for (hsize_t i = 0; i < expect_n; ++i) + rho3d[i] = vp2rho(vp3d[i]); + } + // --- gc / gs (optional, only if model_para_type is MODEL_AZI_ANI) --- if (IP.inversion().model_para_type == MODEL_AZI_ANI) { if (f.exists("gc")) { @@ -462,45 +490,6 @@ void ModelGrid::load_3d_model() { } } } - - // Initialize Vp and density once, after both shear components are ready. - // Radial inversion derives these fields from mean Vs; other modes retain - // the optional input fields and their empirical fallbacks. - if (IP.inversion().model_para_type == MODEL_RADIAL_ANI) { - logger.Info("Computing shared Vp and density from mean Vs for the radial model.", MODULE_GRID); - for (hsize_t i = 0; i < expect_n; ++i) { - vp3d[i] = vs2vp(vsvvsh2vs(vs3d[i], vsh3d[i])); - rho3d[i] = vp2rho(vp3d[i]); - } - } else { - // --- vp (optional, fallback: empirical vs2vp) --- - if (f.exists("vp")) { - interpolate_or_copy("vp", vp3d); - logger.Info( - same_grid ? "Loaded 'vp' from HDF5 file." : - "Loaded and interpolated 'vp' onto current model grid.", - MODULE_GRID - ); - } else { - logger.Info("'vp' not found in HDF5 file, computing from empirical vs2vp.", MODULE_GRID); - for (hsize_t i = 0; i < expect_n; ++i) - vp3d[i] = vs2vp(vs3d[i]); - } - - // --- rho (optional, fallback: empirical vp2rho) --- - if (f.exists("rho")) { - interpolate_or_copy("rho", rho3d); - logger.Info( - same_grid ? "Loaded 'rho' from HDF5 file." : - "Loaded and interpolated 'rho' onto current model grid.", - MODULE_GRID - ); - } else { - logger.Info("'rho' not found in HDF5 file, computing from empirical vp2rho.", MODULE_GRID); - for (hsize_t i = 0; i < expect_n; ++i) - rho3d[i] = vp2rho(vp3d[i]); - } - } } catch (const std::exception &e) { logger.Error(fmt::format( "ModelGrid: failed to load 3D model from HDF5 file '{}': {}", @@ -749,10 +738,10 @@ void ModelGrid::add_radial_aniso_perturbation( // vsh = vsv * sqrt(zeta) auto [vsv_new, vsh_new] = recover_anisotropy(Vs_new, zeta_new); - // vs3d stores Vsv; Vp and density use the mean Vs. + // Update vs3d to store Vs (for consistency with load_3d_model) vs3d[idx] = vsv_new; vsh3d[idx] = vsh_new; - vp3d[idx] = vs2vp(Vs_new); + vp3d[idx] = vs2vp(vsv_new); rho3d[idx] = vp2rho(vp3d[idx]); } } @@ -827,7 +816,7 @@ void ModelGrid::write(const std::string &subname) { f.write_volume("vsv", vs3d, ngrid_i, ngrid_j, ngrid_k); f.write_volume("vsh", vsh3d, ngrid_i, ngrid_j, ngrid_k); - // Compute and write Vs = sqrt((2*vsv^2 + vsh^2)/3) and zeta = vsh^2/vsv^2 + // Compute and write Vs = sqrt((2*vsv + vsh)/3) and zeta = vsh^2/vsv^2 std::vector Vs(nelem); std::vector zeta(nelem); for (int i = 0; i < nelem; ++i) { diff --git a/src/postproc.cpp b/src/postproc.cpp index 147b2ea..b25d8a3 100644 --- a/src/postproc.cpp +++ b/src/postproc.cpp @@ -518,7 +518,7 @@ void postproc::kernel_precondition(SurfGrid& sg) { ); // Precondition the kernels by multiplying with the reference model parameters at each surface grid point - int nker = NPARAMS; + int nker = NPARAMS - 1; for (int iparam = 0; iparam < nker; ++iparam) { if (sg.is_active_ker(iparam)) { sg.ker_loc[iparam] = sg.ker_loc[iparam] * hess_inv; diff --git a/src/preproc.cpp b/src/preproc.cpp index 929696c..5a3fdad 100644 --- a/src/preproc.cpp +++ b/src/preproc.cpp @@ -212,10 +212,6 @@ void preproc::combine_kernels(SurfGrid& sg) { // vs kernel — always allocated (index 0) sg.ker_loc[0] = Tensor3r(dcp.loc_nx(), dcp.loc_ny(), ngrid_k); sg.ker_loc[0].setZero(); - if (IP.inversion().model_para_type == MODEL_RADIAL_ANI) { - sg.ker_loc[5] = Tensor3r(dcp.loc_nx(), dcp.loc_ny(), ngrid_k); - sg.ker_loc[5].setZero(); - } if (IP.inversion().use_alpha_beta_rho && sg.wave_type() == WaveType::RL) { // vp (1) and rho (2) — only when parametrised independently sg.ker_loc[1] = Tensor3r(dcp.loc_nx(), dcp.loc_ny(), ngrid_k); @@ -255,31 +251,7 @@ void preproc::combine_kernels(SurfGrid& sg) { } // Isotropic parameter kernels — gated by use_alpha_beta_rho only - if (IP.inversion().model_para_type == MODEL_RADIAL_ANI) { - for (int ix = 0; ix < dcp.loc_nx(); ++ix) { - for (int iy = 0; iy < dcp.loc_ny(); ++iy) { - const int gx = dcp.loc_I_start() + ix; - const int gy = dcp.loc_J_start() + iy; - for (int k = 0; k < ngrid_k; ++k) { - const real_t vsv = mg.vs3d_loc(ix, iy, k); - const real_t vsh = mg.vsh3d_loc(ix, iy, k); - const auto [ka, kb] = radial_material_kernels( - vsv, vsh, sg.sen_vs_loc(ix, iy, k, iper), - sg.sen_vp_loc(ix, iy, k, iper), - sg.sen_rho_loc(ix, iy, k, iper), - sg.wave_type() == WaveType::LV); - sg.ker_loc[0](ix, iy, k) -= adj_tt(gx, gy) * ka; - sg.ker_loc[5](ix, iy, k) -= adj_tt(gx, gy) * kb; - if (IP.postproc().is_kden) { - // Preserve the wave-specific shear-speed coverage scale, - // including the required empirical-material chain rule. - sg.ker_den_loc(ix, iy, k) -= adj_den(gx, gy) * - (sg.wave_type() == WaveType::LV ? kb : ka); - } - } - } - } - } else if (IP.inversion().use_alpha_beta_rho && sg.wave_type() == WaveType::RL) { + if (IP.inversion().use_alpha_beta_rho && sg.wave_type() == WaveType::RL) { for (int ix = 0; ix < dcp.loc_nx(); ++ix) { for (int iy = 0; iy < dcp.loc_ny(); ++iy) { const int iglob_x = dcp.loc_I_start() + ix; @@ -386,3 +358,4 @@ void preproc::combine_kernels(SurfGrid& sg) { mpi.barrier(); } + diff --git a/src/surf_grid.cpp b/src/surf_grid.cpp index eb3e918..1267bd0 100644 --- a/src/surf_grid.cpp +++ b/src/surf_grid.cpp @@ -82,12 +82,10 @@ SurfGrid::SurfGrid(WaveType wt, SurfType vt){ void SurfGrid::setup_active_kernels() { auto& IP = InputParams::IP(); - active_kernels_ = std::vector(NPARAMS, false); - active_kernels_[0] = true; // vs (Vsv for radial models) - if (IP.inversion().model_para_type == MODEL_RADIAL_ANI) - active_kernels_[5] = true; // independent Vsh partial, for both waves + active_kernels_ = std::vector(NPARAMS-1, false); + active_kernels_[0] = true; // vp if (IP.inversion().use_alpha_beta_rho && wt_ == WaveType::RL) { - active_kernels_[1] = true; // vp + active_kernels_[1] = true; // vs active_kernels_[2] = true; // rho } if (IP.inversion().model_para_type == MODEL_AZI_ANI && wt_ == WaveType::RL) { @@ -300,10 +298,7 @@ void SurfGrid::compute_dispersion_kernel() { vp1d(k) = mg.vp3d_loc(ix, iy, k); rho1d(k) = mg.rho3d_loc(ix, iy, k); } else { - const real_t mean_vs = IP.inversion().model_para_type == MODEL_RADIAL_ANI - ? vsvvsh2vs(mg.vs3d_loc(ix, iy, k), mg.vsh3d_loc(ix, iy, k)) - : vs1d(k); - vp1d(k) = vs2vp(mean_vs); + vp1d(k) = vs2vp(vs1d(k)); rho1d(k) = vp2rho(vp1d(k)); } } @@ -322,11 +317,10 @@ void SurfGrid::compute_dispersion_kernel() { vs_block = kernels.sen_vs.transpose(); if (wt_ == WaveType::RL){ Eigen::Map vp_block(sen_vp_loc.data() + id0, ngrid_k, nperiod_); + Eigen::Map rho_block(sen_rho_loc.data() + id0, ngrid_k, nperiod_); vp_block = kernels.sen_vp.transpose(); + rho_block = kernels.sen_rho.transpose(); } - // Love also has a density kernel; radial empirical scaling needs it. - Eigen::Map rho_block(sen_rho_loc.data() + id0, ngrid_k, nperiod_); - rho_block = kernels.sen_rho.transpose(); if (IP.inversion().model_para_type == MODEL_AZI_ANI) { Eigen::Map gc_block(sen_gc_loc.data() + id0, ngrid_k, nperiod_); Eigen::Map gs_block(sen_gs_loc.data() + id0, ngrid_k, nperiod_); diff --git a/src/surfker/rlker.cpp b/src/surfker/rlker.cpp index acb6d8a..02381f5 100644 --- a/src/surfker/rlker.cpp +++ b/src/surfker/rlker.cpp @@ -790,7 +790,7 @@ static void energy_func(double om, double wvno, RayleighWork& s) double facav = rho * av * av * DUZDUZ / wvno2; s.dcda[m] = facah + facav; double facr = -0.5 * c * c * (URUR + UZUZ); - s.dcdr[m] = 0.5 * (facav + facah) + facr * rho; + s.dcdr[m] = 0.5 * (av * facav + ah * facah) + facr * rho; s.dcdb[m] = 0.0; s.dcdgc[m] = 0.0; s.dcdgs[m] = 0.0; @@ -819,7 +819,7 @@ static void energy_func(double om, double wvno, RayleighWork& s) s.dcda[m] = facah + facav; s.dcdb[m] = facbv + facbh; double facr = -0.5 * c * c * (URUR + UZUZ); - s.dcdr[m] = 0.5 * (facav + facah + facbv) + facr * rho; + s.dcdr[m] = 0.5 * (av * facav + ah * facah + bv * facbv) + facr * rho; s.dcdgc[m] = rho*bv*bv*(UZUZ + 2.0*UZDUR/wvno + DURDUR/wvno2) + rho*bv*bv*URUR; s.dcdgs[m] = s.dcdgc[m]; } @@ -833,12 +833,9 @@ static void energy_func(double om, double wvno, RayleighWork& s) double inv_ug_s0 = 1.0 / (s.ugr * s.sumi0); for (int m = 0; m < mmax; m++) { - // facah/facav/facbv contain rho * velocity^2: dividing only by - // U*I0 gives logarithmic velocity derivatives. The public kernels - // are absolute derivatives dc/dVp, dc/dVs and dc/drho. - s.dcda[m] *= inv_ug_s0 / s.za[m]; - s.dcdb[m] *= (s.iwat[m] == 1) ? 0.0 : inv_ug_s0 / s.zb[m]; - s.dcdr[m] *= inv_ug_s0 / s.zrho[m]; + s.dcda[m] *= inv_ug_s0; + s.dcdb[m] *= inv_ug_s0; + s.dcdr[m] *= inv_ug_s0; s.dcdgc[m] *= inv_ug_s0; s.dcdgs[m] *= inv_ug_s0; } From ea24c086764740839ba12d0bdd22e6fea028e834 Mon Sep 17 00:00:00 2001 From: xumi1993 Date: Sat, 12 Sep 2026 09:47:32 +0800 Subject: [PATCH 4/6] Refactor kernel combination logic for radial anisotropy handling and simplify energy function calculations --- src/optimize.cpp | 13 +++++++------ src/preproc.cpp | 22 +++++++++++++++------- src/surfker/rlker.cpp | 13 ++++++++----- 3 files changed, 30 insertions(+), 18 deletions(-) diff --git a/src/optimize.cpp b/src/optimize.cpp index e346da4..8df13a0 100644 --- a/src/optimize.cpp +++ b/src/optimize.cpp @@ -63,13 +63,14 @@ real_t calc_descent_angle(const FieldVec &direction, const FieldVec &gradient) { // field_dot_global — MPI-reduced directional derivative magnitude // --------------------------------------------------------------------------- // The adjoint kernels are densities over the horizontal surface, while their -// depth dimension already consists of discrete layer sensitivities (dc/dm_k). -// The search direction, however, is applied multiplicatively to vs/vp/rho/gamma -// and additively to gc/gs. Therefore this routine includes both the horizontal -// quadrature weight and the chain rule from alpha to the physical model: +// depth dimension already consists of discrete layer sensitivities. This +// routine therefore adds only the horizontal quadrature weight. // -// m(alpha) = m0 (1 - alpha*d) => -dm/dalpha = m0*d -// m(alpha) = m0 - alpha*d => -dm/dalpha = d . +// No per-parameter chain rule from alpha to the physical model is needed: +// combine_kernels() already returns dchi/dln(m) for the multiplicatively +// updated parameters (vs/vp/rho/gamma) and dchi/dm for the additive ones +// (gc/gs), so in both cases the step is x(alpha) = x0 - alpha*d and +// -dx/dalpha = d. // // dlon_i and dlat_j are nodal quadrature widths. End points receive half of // the adjacent interval, as in the trapezoidal rule. No dz factor is added. diff --git a/src/preproc.cpp b/src/preproc.cpp index 5a3fdad..4f93962 100644 --- a/src/preproc.cpp +++ b/src/preproc.cpp @@ -289,17 +289,25 @@ void preproc::combine_kernels(SurfGrid& sg) { const real_t vp = mg.vp3d_loc(ix, iy, k); const real_t dab = dalpha_dbeta(vs); // d(vp)/d(vs) const real_t dra = drho_dalpha(vp); // d(rho)/d(vp) - sg.ker_loc[0](ix, iy, k) -= att * ( - sg.sen_vs_loc(ix, iy, k, iper) - + sg.sen_vp_loc(ix, iy, k, iper) * dab - + sg.sen_rho_loc(ix, iy, k, iper) * dra * dab - ); - if (IP.postproc().is_kden) { - sg.ker_den_loc(ix, iy, k) -= adj_den(iglob_x, iglob_y) * ( + if ( IP.inversion().model_para_type == MODEL_RADIAL_ANI) { + sg.ker_loc[0](ix, iy, k) -= att * sg.sen_vs_loc(ix, iy, k, iper); + if (IP.postproc().is_kden) { + sg.ker_den_loc(ix, iy, k) -= adj_den(iglob_x, iglob_y) * sg.sen_vs_loc(ix, iy, k, iper); + } + } else { + sg.ker_loc[0](ix, iy, k) -= att * ( sg.sen_vs_loc(ix, iy, k, iper) + sg.sen_vp_loc(ix, iy, k, iper) * dab + sg.sen_rho_loc(ix, iy, k, iper) * dra * dab ); + + if (IP.postproc().is_kden) { + sg.ker_den_loc(ix, iy, k) -= adj_den(iglob_x, iglob_y) * ( + sg.sen_vs_loc(ix, iy, k, iper) + + sg.sen_vp_loc(ix, iy, k, iper) * dab + + sg.sen_rho_loc(ix, iy, k, iper) * dra * dab + ); + } } } } diff --git a/src/surfker/rlker.cpp b/src/surfker/rlker.cpp index 02381f5..acb6d8a 100644 --- a/src/surfker/rlker.cpp +++ b/src/surfker/rlker.cpp @@ -790,7 +790,7 @@ static void energy_func(double om, double wvno, RayleighWork& s) double facav = rho * av * av * DUZDUZ / wvno2; s.dcda[m] = facah + facav; double facr = -0.5 * c * c * (URUR + UZUZ); - s.dcdr[m] = 0.5 * (av * facav + ah * facah) + facr * rho; + s.dcdr[m] = 0.5 * (facav + facah) + facr * rho; s.dcdb[m] = 0.0; s.dcdgc[m] = 0.0; s.dcdgs[m] = 0.0; @@ -819,7 +819,7 @@ static void energy_func(double om, double wvno, RayleighWork& s) s.dcda[m] = facah + facav; s.dcdb[m] = facbv + facbh; double facr = -0.5 * c * c * (URUR + UZUZ); - s.dcdr[m] = 0.5 * (av * facav + ah * facah + bv * facbv) + facr * rho; + s.dcdr[m] = 0.5 * (facav + facah + facbv) + facr * rho; s.dcdgc[m] = rho*bv*bv*(UZUZ + 2.0*UZDUR/wvno + DURDUR/wvno2) + rho*bv*bv*URUR; s.dcdgs[m] = s.dcdgc[m]; } @@ -833,9 +833,12 @@ static void energy_func(double om, double wvno, RayleighWork& s) double inv_ug_s0 = 1.0 / (s.ugr * s.sumi0); for (int m = 0; m < mmax; m++) { - s.dcda[m] *= inv_ug_s0; - s.dcdb[m] *= inv_ug_s0; - s.dcdr[m] *= inv_ug_s0; + // facah/facav/facbv contain rho * velocity^2: dividing only by + // U*I0 gives logarithmic velocity derivatives. The public kernels + // are absolute derivatives dc/dVp, dc/dVs and dc/drho. + s.dcda[m] *= inv_ug_s0 / s.za[m]; + s.dcdb[m] *= (s.iwat[m] == 1) ? 0.0 : inv_ug_s0 / s.zb[m]; + s.dcdr[m] *= inv_ug_s0 / s.zrho[m]; s.dcdgc[m] *= inv_ug_s0; s.dcdgs[m] *= inv_ug_s0; } From 85794d55630dbbd94db4e123e79661f7ece17bc6 Mon Sep 17 00:00:00 2001 From: Mijian Xu Date: Mon, 14 Sep 2026 15:57:29 +0800 Subject: [PATCH 5/6] Refactor kernel adjustment logic for anisotropic model Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/preproc.cpp | 23 ++++++++--------------- 1 file changed, 8 insertions(+), 15 deletions(-) diff --git a/src/preproc.cpp b/src/preproc.cpp index 4f93962..40e5d07 100644 --- a/src/preproc.cpp +++ b/src/preproc.cpp @@ -289,25 +289,18 @@ void preproc::combine_kernels(SurfGrid& sg) { const real_t vp = mg.vp3d_loc(ix, iy, k); const real_t dab = dalpha_dbeta(vs); // d(vp)/d(vs) const real_t dra = drho_dalpha(vp); // d(rho)/d(vp) - if ( IP.inversion().model_para_type == MODEL_RADIAL_ANI) { - sg.ker_loc[0](ix, iy, k) -= att * sg.sen_vs_loc(ix, iy, k, iper); - if (IP.postproc().is_kden) { - sg.ker_den_loc(ix, iy, k) -= adj_den(iglob_x, iglob_y) * sg.sen_vs_loc(ix, iy, k, iper); - } - } else { - sg.ker_loc[0](ix, iy, k) -= att * ( + sg.ker_loc[0](ix, iy, k) -= att * ( + sg.sen_vs_loc(ix, iy, k, iper) + + sg.sen_vp_loc(ix, iy, k, iper) * dab + + sg.sen_rho_loc(ix, iy, k, iper) * dra * dab + ); + + if (IP.postproc().is_kden) { + sg.ker_den_loc(ix, iy, k) -= adj_den(iglob_x, iglob_y) * ( sg.sen_vs_loc(ix, iy, k, iper) + sg.sen_vp_loc(ix, iy, k, iper) * dab + sg.sen_rho_loc(ix, iy, k, iper) * dra * dab ); - - if (IP.postproc().is_kden) { - sg.ker_den_loc(ix, iy, k) -= adj_den(iglob_x, iglob_y) * ( - sg.sen_vs_loc(ix, iy, k, iper) - + sg.sen_vp_loc(ix, iy, k, iper) * dab - + sg.sen_rho_loc(ix, iy, k, iper) * dra * dab - ); - } } } } From fb02330f7e6a6734386b3423641fa1a04ec24f3c Mon Sep 17 00:00:00 2001 From: Mijian Xu Date: Mon, 14 Sep 2026 15:57:41 +0800 Subject: [PATCH 6/6] Refactor comments in optimize.cpp for clarity Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com> --- src/optimize.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/optimize.cpp b/src/optimize.cpp index 8df13a0..c1fc141 100644 --- a/src/optimize.cpp +++ b/src/optimize.cpp @@ -66,11 +66,10 @@ real_t calc_descent_angle(const FieldVec &direction, const FieldVec &gradient) { // depth dimension already consists of discrete layer sensitivities. This // routine therefore adds only the horizontal quadrature weight. // -// No per-parameter chain rule from alpha to the physical model is needed: -// combine_kernels() already returns dchi/dln(m) for the multiplicatively -// updated parameters (vs/vp/rho/gamma) and dchi/dm for the additive ones -// (gc/gs), so in both cases the step is x(alpha) = x0 - alpha*d and -// -dx/dalpha = d. +// The line-search derivative includes the model-update chain rule below: +// for multiplicative parameters, -dm/dalpha = m*d; for additive gc/gs, +// -dm/dalpha = d. The directional derivative therefore uses the physical +// model values for vs/vp/rho and gamma before applying the horizontal weight. // // dlon_i and dlat_j are nodal quadrature widths. End points receive half of // the adjacent interval, as in the trapezoidal rule. No dz factor is added.