168 lines
6.5 KiB
C++
168 lines
6.5 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* 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 <http://www.gnu.org/licenses/>.
|
|
*
|
|
* 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 "mkernel_block.h"
|
|
|
|
#ifdef GCTL_POTENTIAL_AUTODIFF
|
|
|
|
autodiff::var magobser_block_thickness_za_sig(const gctl::array<gctl::mag_block> &ele, const gctl::array<double> &sus,
|
|
const gctl::point3dc &obsp, const autodiff::ArrayXvar& btm);
|
|
|
|
double gctl::magobser_wrt_thickness(array<double> &btm_derivatives, loss_func &lf,
|
|
const array<point3dc> &obsp, const array<double> &sus, const array<mag_block> &ele,
|
|
magnetic_field_type_e comp_id, verbose_type_e verbose)
|
|
{
|
|
if (ele[0].att == nullptr) throw gctl::runtime_error("[GCTL::Magkernel_Block] Magnetization parameter not set.");
|
|
|
|
int i, j;
|
|
int o_size = obsp.size();
|
|
int e_size = ele.size();
|
|
|
|
btm_derivatives.resize(e_size, 0.0);
|
|
|
|
double tmp_val, tmp_grad;
|
|
autodiff::var tmp_obs;
|
|
autodiff::ArrayXvar btm_depth;
|
|
Eigen::VectorXd btm_ders;
|
|
|
|
gctl::progress_bar bar(o_size, "Za_wrt_thickness (openmp)");
|
|
if (verbose == gctl::FullMsg) bar.progressed(0);
|
|
else if (verbose == gctl::ShortMsg) bar.progressed_simple(0);
|
|
|
|
#pragma omp parallel for private (i, j, tmp_val, tmp_grad, tmp_obs, btm_depth, btm_ders) shared(e_size) schedule(guided)
|
|
for (i = 0; i < o_size; i++)
|
|
{
|
|
btm_depth.resize(e_size);
|
|
for (j = 0; j < e_size; j++)
|
|
{
|
|
btm_depth[j] = ele[j].dl->z;
|
|
}
|
|
|
|
btm_ders.setZero();
|
|
tmp_obs = magobser_block_thickness_za_sig(ele, sus, obsp[i], btm_depth);
|
|
btm_ders = autodiff::gradient(tmp_obs, btm_depth);
|
|
tmp_val = autodiff::val(tmp_obs);
|
|
tmp_grad = lf.gradient(tmp_val, i);
|
|
|
|
#pragma omp critical
|
|
{
|
|
lf.evaluate(tmp_val, i);
|
|
for (j = 0; j < e_size; j++)
|
|
{
|
|
btm_derivatives[j] += tmp_grad*btm_ders[j];
|
|
}
|
|
}
|
|
//if (verbose == gctl::FullMsg) bar.progressed(i);
|
|
//else if (verbose == gctl::ShortMsg) bar.progressed_simple(i);
|
|
}
|
|
|
|
if (verbose == gctl::FullMsg) bar.progressed(o_size - 1);
|
|
else if (verbose == gctl::ShortMsg) bar.progressed_simple(o_size - 1);
|
|
|
|
return lf.get_loss();
|
|
}
|
|
|
|
autodiff::var magobser_block_thickness_za_sig(const gctl::array<gctl::mag_block> &ele, const gctl::array<double> &sus,
|
|
const gctl::point3dc &obsp, const autodiff::ArrayXvar& btm)
|
|
{
|
|
double I,A;
|
|
double Alpha,Beta,Gamma;
|
|
double x, y, z, su;
|
|
double x1,x2,y1,y2,z2;
|
|
double w111_1,w112_1,w121_1,w122_1,w211_1,w212_1,w221_1,w222_1;
|
|
double w222_2,w112_2,w122_2,w212_2;
|
|
double R222,R122,R212,R112;
|
|
double G222,G122,G212,G112;
|
|
|
|
autodiff::var obs = 0.0;
|
|
autodiff::var w211_2,w121_2,w221_2,w111_2;
|
|
autodiff::var R221,R121,R211,R111;
|
|
autodiff::var G221,G121,G211,G111;
|
|
|
|
gctl::magblock_para* mp;
|
|
for (int i = 0; i < ele.size(); i++)
|
|
{
|
|
mp = ele[i].att;
|
|
I = mp->inclina_deg*M_PI/180.0;
|
|
A = (90 - mp->declina_deg)*M_PI/180.0;
|
|
|
|
Alpha=cos(I)*cos(A);
|
|
Beta=cos(I)*sin(A);
|
|
Gamma=sin(I);
|
|
|
|
x = obsp.x;
|
|
y = obsp.y;
|
|
z = obsp.z;
|
|
su = sus[i];
|
|
|
|
x1 = ele[i].dl->x; x2 = ele[i].ur->x;
|
|
y1 = ele[i].dl->y; y2 = ele[i].ur->y;
|
|
z2 = ele[i].ur->z;
|
|
|
|
R222=sqrt((x2-x)*(x2-x)+(y2-y)*(y2-y)+(z2-z)*(z2-z));
|
|
R122=sqrt((x1-x)*(x1-x)+(y2-y)*(y2-y)+(z2-z)*(z2-z));
|
|
R212=sqrt((x2-x)*(x2-x)+(y1-y)*(y1-y)+(z2-z)*(z2-z));
|
|
R112=sqrt((x1-x)*(x1-x)+(y1-y)*(y1-y)+(z2-z)*(z2-z));
|
|
R221=sqrt((x2-x)*(x2-x)+(y2-y)*(y2-y)+(btm[i] - z)*(btm[i] - z));
|
|
R121=sqrt((x1-x)*(x1-x)+(y2-y)*(y2-y)+(btm[i] - z)*(btm[i] - z));
|
|
R211=sqrt((x2-x)*(x2-x)+(y1-y)*(y1-y)+(btm[i] - z)*(btm[i] - z));
|
|
R111=sqrt((x1-x)*(x1-x)+(y1-y)*(y1-y)+(btm[i] - z)*(btm[i] - z));
|
|
|
|
w222_1=(x2-x)*(y2-y);
|
|
w122_1=(x1-x)*(y2-y);
|
|
w212_1=(x2-x)*(y1-y);
|
|
w112_1=(x1-x)*(y1-y);
|
|
w221_1=(x2-x)*(y2-y);
|
|
w121_1=(x1-x)*(y2-y);
|
|
w211_1=(x2-x)*(y1-y);
|
|
w111_1=(x1-x)*(y1-y);
|
|
|
|
w222_2=R222*(z2-z);
|
|
w122_2=R122*(z2-z);
|
|
w212_2=R212*(z2-z);
|
|
w112_2=R112*(z2-z);
|
|
w221_2=R221*(btm[i] - z);
|
|
w121_2=R121*(btm[i] - z);
|
|
w211_2=R211*(btm[i] - z);
|
|
w111_2=R111*(btm[i] - z);
|
|
|
|
G222=-Gamma*atan2(w222_1,w222_2)+Alpha*log(R222+y2-y)+Beta*log(R222+x2-x);
|
|
G122=-Gamma*atan2(w122_1,w122_2)+Alpha*log(R122+y2-y)+Beta*log(R122+x1-x);
|
|
G212=-Gamma*atan2(w212_1,w212_2)+Alpha*log(R212+y1-y)+Beta*log(R212+x2-x);
|
|
G112=-Gamma*atan2(w112_1,w112_2)+Alpha*log(R112+y1-y)+Beta*log(R112+x1-x);
|
|
G221=-Gamma*atan2(w221_1,w221_2)+Alpha*log(R221+y2-y)+Beta*log(R221+x2-x);
|
|
G121=-Gamma*atan2(w121_1,w121_2)+Alpha*log(R121+y2-y)+Beta*log(R121+x1-x);
|
|
G211=-Gamma*atan2(w211_1,w211_2)+Alpha*log(R211+y1-y)+Beta*log(R211+x2-x);
|
|
G111=-Gamma*atan2(w111_1,w111_2)+Alpha*log(R111+y1-y)+Beta*log(R111+x1-x);
|
|
|
|
obs += GCTL_T0*su/(4*M_PI)*(G222-G122-G212+G112-G221+G121+G211-G111);
|
|
}
|
|
|
|
return obs;
|
|
}
|
|
|
|
#endif // GCTL_POTENTIAL_AUTODIFF
|