Devel

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)
 

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();

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<int>(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();

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.

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)

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;
     }

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 << "      </Grid>\n";
     }