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/optimize.cpp b/src/optimize.cpp index e346da4..c1fc141 100644 --- a/src/optimize.cpp +++ b/src/optimize.cpp @@ -63,13 +63,13 @@ 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 . +// 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. diff --git a/src/preproc.cpp b/src/preproc.cpp index 5a3fdad..40e5d07 100644 --- a/src/preproc.cpp +++ b/src/preproc.cpp @@ -294,6 +294,7 @@ void preproc::combine_kernels(SurfGrid& sg) { + 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) 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; } 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"; }