diff --git a/CMakeLists.txt b/CMakeLists.txt index ea09cec..0d6781c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -84,5 +84,5 @@ configure_file( # 添加库源文件地址 add_subdirectory(lib) -add_subdirectory(tool) add_subdirectory(example) +add_subdirectory(tool) \ No newline at end of file diff --git a/lib/potential/gkernel_sphere.h b/lib/potential/gkernel_sphere.h index 57495d0..6f50be6 100644 --- a/lib/potential/gkernel_sphere.h +++ b/lib/potential/gkernel_sphere.h @@ -31,9 +31,15 @@ #include "gctl/core/array.h" #include "gctl/geometry/point3c.h" #include "gctl/geometry/sphere.h" +#include "gctl/optimization/loss_func.h" #include "gm_data.h" +#ifdef GCTL_POTENTIAL_AUTODIFF +#include "autodiff/reverse/var.hpp" +#include "autodiff/reverse/var/eigen.hpp" +#endif // GCTL_POTENTIAL_AUTODIFF + namespace gctl { /** @@ -60,6 +66,24 @@ namespace gctl */ void gobser_sphere(array &out_obs, const array &obsp, const sphere &gs, gravitational_field_type_e comp_id = Vz, value_operator_e valtype = ReplaceVal); + +#ifdef GCTL_POTENTIAL_AUTODIFF + + /** + * @brief 计算球体的正演重力异常与观测数据的数据拟合差,以及其相对于球体质量与埋深的偏导数 + * + * @param depth_der 相对于埋深的偏导数 + * @param mass_der 相对于质量的偏导数 + * @param mass_loc 球体的位置 + * @param mass 球体的质量 + * @param lf 数据拟合差对象(包含了观测数据与不确定度等信息) + * @param obsp 观测点位 + */ + void gobser_wrt_mass_depth(double &depth_der, double &mass_der, const point3dc &mass_loc, + double mass, loss_func &lf, const array &obsp); + +#endif // GCTL_POTENTIAL_AUTODIFF + } #endif // _GCTL_GRAV_KERNEL_SPHERE_H \ No newline at end of file diff --git a/lib/potential/gkernel_sphere_mass_depth.cpp b/lib/potential/gkernel_sphere_mass_depth.cpp new file mode 100644 index 0000000..8b6dc92 --- /dev/null +++ b/lib/potential/gkernel_sphere_mass_depth.cpp @@ -0,0 +1,83 @@ +/******************************************************** + * ██████╗ ██████╗████████╗██╗ + * ██╔════╝ ██╔════╝╚══██╔══╝██║ + * ██║ ███╗██║ ██║ ██║ + * ██║ ██║██║ ██║ ██║ + * ╚██████╔╝╚██████╗ ██║ ███████╗ + * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ + * Geophysical Computational Tools & Library (GCTL) + * + * Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn) + * + * GCTL is distributed under a dual licensing scheme. You can redistribute + * it and/or modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation, either version 2 + * of the License, or (at your option) any later version. You should have + * received a copy of the GNU Lesser General Public License along with this + * program. If not, see . + * + * If the terms and conditions of the LGPL v.2. would prevent you from using + * the GCTL, please consider the option to obtain a commercial license for a + * fee. These licenses are offered by the GCTL's original author. As a rule, + * licenses are provided "as-is", unlimited in time for a one time fee. Please + * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget + * to include some description of your company and the realm of its activities. + * Also add information on how to contact you by electronic and paper mail. + ******************************************************/ + +#include "gkernel_sphere.h" + +#ifdef GCTL_POTENTIAL_AUTODIFF + +autodiff::var gobser_mass_depth_vz_sig(const gctl::point3dc &obsp, double loc_x, + double loc_y, const autodiff::ArrayXvar& depth_mass); + +void gctl::gobser_wrt_mass_depth(double &depth_der, double &mass_der, const point3dc &mass_loc, + double mass, loss_func &lf, const array &obsp) +{ + double tmp_val, tmp_grad; + autodiff::var tmp_obs; + autodiff::ArrayXvar depth_mass; + Eigen::VectorXd ders; + + depth_der = 0.0; + mass_der = 0.0; + + int i; + int o_size = obsp.size(); + +#pragma omp parallel for private (i, tmp_val, tmp_grad, tmp_obs, depth_mass, ders) schedule(guided) + for (i = 0; i < o_size; i++) + { + depth_mass.resize(2); + depth_mass[0] = mass_loc.z; + depth_mass[1] = mass; + + ders.setZero(); + tmp_obs = gobser_mass_depth_vz_sig(obsp[i], mass_loc.x, mass_loc.y, depth_mass); + ders = autodiff::gradient(tmp_obs, depth_mass); + tmp_val = autodiff::val(tmp_obs); + tmp_grad = lf.gradient(tmp_val, i); + +#pragma omp critical +{ + lf.evaluate(tmp_val, i); + depth_der += tmp_grad*ders[0]; + mass_der += tmp_grad*ders[1]; +} + } + + return; +} + + +autodiff::var gobser_mass_depth_vz_sig(const gctl::point3dc &obsp, double loc_x, + double loc_y, const autodiff::ArrayXvar& depth_mass) +{ + autodiff::var tmp_obs = 1e+8*GCTL_G0*depth_mass[1]*(obsp.z - depth_mass[0]) + /pow(pow(obsp.x - loc_x, 2) + pow(obsp.y - loc_y, 2) + pow(obsp.z - depth_mass[0], 2), 1.5); + + return tmp_obs; +} + +#endif // GCTL_POTENTIAL_AUTODIFF \ No newline at end of file diff --git a/lib/potential/mkernel_block.h b/lib/potential/mkernel_block.h index 3e1d3fc..a7afb35 100644 --- a/lib/potential/mkernel_block.h +++ b/lib/potential/mkernel_block.h @@ -35,6 +35,7 @@ #include "autodiff/reverse/var.hpp" #include "autodiff/reverse/var/eigen.hpp" #endif // GCTL_POTENTIAL_AUTODIFF + namespace gctl { struct magblock_para