From 4db4588993ded6b241eeee68b988e281afe3a4d0 Mon Sep 17 00:00:00 2001
From: XavierAdriaens <xavier.adriaens@doct.uliege.be>
Date: Fri, 9 Dec 2022 23:15:10 +0100
Subject: [PATCH] Phase only objective function
---
common/data/element.cpp | 36 +++++++++++++++++
common/data/element.h | 25 ++++++++++--
specific/data/element.cpp | 30 ++++++++++++++
specific/data/element.h | 8 ++++
specific/data/objective/newObjective.h | 2 +
specific/data/objective/phase.cpp | 54 +++++++++++++++++++++++++
specific/data/objective/phase.h | 29 +++++++++++++
specific/model/innerproduct/sobolev.cpp | 14 +++----
8 files changed, 188 insertions(+), 10 deletions(-)
create mode 100644 specific/data/objective/phase.cpp
create mode 100644 specific/data/objective/phase.h
diff --git a/common/data/element.cpp b/common/data/element.cpp
index 7195ef0..3f281fc 100644
--- a/common/data/element.cpp
+++ b/common/data/element.cpp
@@ -266,6 +266,42 @@ Data<T_Physic> Data<T_Physic>::operator-(const Data<T_Physic>& other) const
return sub;
}
+template<Physic T_Physic>
+Data<T_Physic> Data<T_Physic>::operator*(const Data<T_Physic>& other) const
+{
+ areCompatible(*this,other);
+ Data<T_Physic> out(*this);
+ for (unsigned int f = 0; f < nf(); f++)
+ {
+ for (unsigned int s = 0; s < ns(); s++)
+ {
+ for (unsigned int r = 0; r < nr(s); r++)
+ {
+ out.value(f,s,r, wise_multiply<T_Physic>(this->value(f,s,r),other.value(f,s,r)) );
+ }
+ }
+ }
+ return out;
+}
+
+template<Physic T_Physic>
+Data<T_Physic> Data<T_Physic>::operator/(const Data<T_Physic>& other) const
+{
+ areCompatible(*this,other);
+ Data<T_Physic> out(*this);
+ for (unsigned int f = 0; f < nf(); f++)
+ {
+ for (unsigned int s = 0; s < ns(); s++)
+ {
+ for (unsigned int r = 0; r < nr(s); r++)
+ {
+ out.value(f,s,r, wise_divide<T_Physic>(this->value(f,s,r),other.value(f,s,r)) );
+ }
+ }
+ }
+ return out;
+}
+
template<Physic T_Physic>
Data<T_Physic> Data<T_Physic>::operator*(std::complex<double> a) const
{
diff --git a/common/data/element.h b/common/data/element.h
index cd3667f..1198126 100644
--- a/common/data/element.h
+++ b/common/data/element.h
@@ -43,6 +43,8 @@ public:
void operator=(const Data<T_Physic>& other);
Data<T_Physic> operator+(const Data<T_Physic>& other)const;
Data<T_Physic> operator-(const Data<T_Physic>& other)const;
+ Data<T_Physic> operator*(const Data<T_Physic>& other)const;
+ Data<T_Physic> operator/(const Data<T_Physic>& other)const;
Data<T_Physic> operator*(std::complex<double> a)const;
Data<T_Physic> operator/(std::complex<double> a)const;
private:
@@ -132,18 +134,35 @@ Data<T_Physic> conj(const Data<T_Physic>& d)
template<Physic T_Physic>
Data<T_Physic> abs(const Data<T_Physic>& d)
{
- Data<T_Physic> abs(d);
+ Data<T_Physic> out(d);
for (unsigned int f = 0; f < d.nf(); f++)
{
for (unsigned int s = 0; s < d.ns(); s++)
{
for (unsigned int r = 0; r < d.nr(s); r++)
{
- abs.value(f,s,r, std::complex<double>(std::abs(d.value(f,s,r)) ) );
+ out.value(f,s,r, abs<T_Physic>(d.value(f,s,r)) );
}
}
}
- return abs;
+ return out;
+};
+
+template<Physic T_Physic>
+Data<T_Physic> real(const Data<T_Physic>& d)
+{
+ Data<T_Physic> out(d);
+ for (unsigned int f = 0; f < d.nf(); f++)
+ {
+ for (unsigned int s = 0; s < d.ns(); s++)
+ {
+ for (unsigned int r = 0; r < d.nr(s); r++)
+ {
+ out.value(f,s,r, real<T_Physic>(d.value(f,s,r)) );
+ }
+ }
+ }
+ return out;
};
#endif // H_COMMON_DATA_ELEMENT
diff --git a/specific/data/element.cpp b/specific/data/element.cpp
index 53e0513..03bc76a 100644
--- a/specific/data/element.cpp
+++ b/specific/data/element.cpp
@@ -15,6 +15,14 @@ double norm2<Physic::acoustic>(const typename PointData<Physic::acoustic>::Objec
template<>
typename PointData<Physic::acoustic>::Object conj<Physic::acoustic>(const typename PointData<Physic::acoustic>::Object& d){return std::conj(d);}
template<>
+typename PointData<Physic::acoustic>::Object abs<Physic::acoustic>(const typename PointData<Physic::acoustic>::Object& d){return std::abs(d);}
+template<>
+typename PointData<Physic::acoustic>::Object real<Physic::acoustic>(const typename PointData<Physic::acoustic>::Object& d){return std::real(d);}
+template<>
+typename PointData<Physic::acoustic>::Object wise_divide<Physic::acoustic>(const typename PointData<Physic::acoustic>::Object& d1,const typename PointData<Physic::acoustic>::Object& d2){return d1/d2;}
+template<>
+typename PointData<Physic::acoustic>::Object wise_multiply<Physic::acoustic>(const typename PointData<Physic::acoustic>::Object& d1,const typename PointData<Physic::acoustic>::Object& d2){return d1*d2;}
+template<>
typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::acoustic>::degree >::Object s_ones<Physic::acoustic>(){return 1.;}
/* electromagnetic */
@@ -23,6 +31,17 @@ double norm2<Physic::electromagnetic>(const typename PointData<Physic::electroma
template<>
typename PointData<Physic::electromagnetic>::Object conj<Physic::electromagnetic>(const typename PointData<Physic::electromagnetic>::Object& d){return d.conjugate();}
template<>
+typename PointData<Physic::electromagnetic>::Object real<Physic::electromagnetic>(const typename PointData<Physic::electromagnetic>::Object& d){return d.real();}
+template<>
+typename PointData<Physic::electromagnetic>::Object abs<Physic::electromagnetic>(const typename PointData<Physic::electromagnetic>::Object& d)
+{
+ return PointData<Physic::electromagnetic>::ones() * d.norm();
+}
+template<>
+typename PointData<Physic::electromagnetic>::Object wise_divide<Physic::electromagnetic>(const typename PointData<Physic::electromagnetic>::Object& d1,const typename PointData<Physic::electromagnetic>::Object& d2){return d1.array() / d2.array() ;}
+template<>
+typename PointData<Physic::electromagnetic>::Object wise_multiply<Physic::electromagnetic>(const typename PointData<Physic::electromagnetic>::Object& d1,const typename PointData<Physic::electromagnetic>::Object& d2){return d1.array() * d2.array() ;}
+template<>
typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::electromagnetic>::degree >::Object s_ones<Physic::electromagnetic>(){return PointData<Physic::electromagnetic>::Object(1.,1.,1.);}
/* elastodynamic */
@@ -31,4 +50,15 @@ double norm2<Physic::elastodynamic>(const typename PointData<Physic::elastodynam
template<>
typename PointData<Physic::elastodynamic>::Object conj<Physic::elastodynamic>(const typename PointData<Physic::elastodynamic>::Object& d){return d.conjugate();}
template<>
+typename PointData<Physic::elastodynamic>::Object real<Physic::elastodynamic>(const typename PointData<Physic::elastodynamic>::Object& d){return d.real();}
+template<>
+typename PointData<Physic::elastodynamic>::Object abs<Physic::elastodynamic>(const typename PointData<Physic::elastodynamic>::Object& d)
+{
+ return PointData<Physic::elastodynamic>::ones() * d.norm();
+}
+template<>
+typename PointData<Physic::elastodynamic>::Object wise_divide<Physic::elastodynamic>(const typename PointData<Physic::elastodynamic>::Object& d1,const typename PointData<Physic::elastodynamic>::Object& d2){return d1.array() / d2.array() ;}
+template<>
+typename PointData<Physic::elastodynamic>::Object wise_multiply<Physic::elastodynamic>(const typename PointData<Physic::elastodynamic>::Object& d1,const typename PointData<Physic::elastodynamic>::Object& d2){return d1.array() * d2.array() ;}
+template<>
typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::elastodynamic>::degree >::Object s_ones<Physic::elastodynamic>(){return PointData<Physic::elastodynamic>::Object(1.,1.,1.);}
diff --git a/specific/data/element.h b/specific/data/element.h
index c70de14..45654f8 100644
--- a/specific/data/element.h
+++ b/specific/data/element.h
@@ -28,5 +28,13 @@ template<Physic T_Physic>
double norm2(const typename PointData<T_Physic>::Object& d);
template<Physic T_Physic>
typename PointData<T_Physic>::Object conj(const typename PointData<T_Physic>::Object& d);
+template<Physic T_Physic>
+typename PointData<T_Physic>::Object real(const typename PointData<T_Physic>::Object& d);
+template<Physic T_Physic>
+typename PointData<T_Physic>::Object abs(const typename PointData<T_Physic>::Object& d);
+template<Physic T_Physic>
+typename PointData<T_Physic>::Object wise_divide(const typename PointData<T_Physic>::Object& d1,const typename PointData<T_Physic>::Object& d2);
+template<Physic T_Physic>
+typename PointData<T_Physic>::Object wise_multiply(const typename PointData<T_Physic>::Object& d1,const typename PointData<T_Physic>::Object& d2);
#endif // H_SPECIFIC_DATA_ELEMENT
diff --git a/specific/data/objective/newObjective.h b/specific/data/objective/newObjective.h
index 0d57e8a..cb6a8a8 100644
--- a/specific/data/objective/newObjective.h
+++ b/specific/data/objective/newObjective.h
@@ -9,6 +9,7 @@
//GmshFWI Library
#include "l2distance.h"
+#include "phase.h"
template<Physic T_Physic>
ObjectiveInterface<T_Physic>* newObjective(const Data<T_Physic>& d0, const ConfigurationInterface* const config, const gmshfem::common::GmshFem& gmshFem,std::string suffix = "")
@@ -19,6 +20,7 @@ ObjectiveInterface<T_Physic>* newObjective(const Data<T_Physic>& d0, const Confi
throw gmshfem::common::Exception("Objective type could not be found.");
}
if(objective=="l2distance"){return new l2distance::Objective<T_Physic>(d0,gmshFem);}
+ else if(objective=="phase"){return new phase::Objective<T_Physic>(d0,gmshFem);}
else
{
throw gmshfem::common::Exception("Objective " + objective + " is not valid.");
diff --git a/specific/data/objective/phase.cpp b/specific/data/objective/phase.cpp
new file mode 100644
index 0000000..8df2b56
--- /dev/null
+++ b/specific/data/objective/phase.cpp
@@ -0,0 +1,54 @@
+//GmshFem Library
+#include "Exception.h"
+
+//GmshFEM Library
+#include "phase.h"
+
+#define d0 _d0
+#define d ds.state(Type::FS)
+#define rho ds.state(Type::AS)
+#define dd ds.state(Type::PFS)
+#define drho ds.state(Type::PAS)
+
+using namespace gmshfem::common;
+/*
+* PHASE
+*/
+namespace phase
+{
+ template<Physic T_Physic>
+ double Objective<T_Physic>::performance(const Data<T_Physic>& d1)
+ {
+ Data<T_Physic> r = d1/abs(d1)-d0/abs(d0);
+ return 0.5 * l2norm2(r);
+ }
+
+ template<Physic T_Physic>
+ const Data<T_Physic>& Objective<T_Physic>::update(Type type, const DataStateEvaluator<T_Physic>& ds)
+ {
+ switch (type)
+ {
+ case Type::FS: case Type::PFS: default:
+ throw Exception("Objective can not update (perturbed) forward data.");
+ break;
+ case Type::AS:
+ {
+ Data<T_Physic> absd = abs(d);
+ Data<T_Physic> cnjr = conj(d/absd-d0/abs(d0));
+ _v = cnjr / absd - real(cnjr*d) * conj(d) / (absd*absd*absd);
+ break;
+ }
+ case Type::PAS:
+ throw Exception("Objective perturbed adjoint state for phase is not defined. Linear/Anti-linear problem.");
+ break;
+ case Type::DS:
+ /* Does not take into account user sparsity input */
+ _v.value( PointData<T_Physic>::ones() );
+ break;
+ }
+ return _v;
+ }
+ template class Objective<Physic::acoustic>;
+ template class Objective<Physic::electromagnetic>;
+ template class Objective<Physic::elastodynamic>;
+}
diff --git a/specific/data/objective/phase.h b/specific/data/objective/phase.h
new file mode 100644
index 0000000..c653df7
--- /dev/null
+++ b/specific/data/objective/phase.h
@@ -0,0 +1,29 @@
+#ifndef H_SPECIFIC_DATA_OBJECTIVE_PHASE
+#define H_SPECIFIC_DATA_OBJECTIVE_PHASE
+
+//GmshFEM Library
+#include "GmshFem.h"
+
+//GmshFWI Library
+#include "../../../common/data/objective/objective.h"
+
+/*
+* PHASE
+*/
+namespace phase
+{
+ template<Physic T_Physic>
+ class Objective final: public ObjectiveInterface<T_Physic>
+ {
+ private:
+ using ObjectiveInterface<T_Physic>::_d0;
+ using ObjectiveInterface<T_Physic>::_v;
+ public:
+ Objective(const Data<T_Physic>& d0, const gmshfem::common::GmshFem& gmshFem) : ObjectiveInterface<T_Physic>(d0) {};
+
+ virtual double performance(const Data<T_Physic>& d);
+ virtual const Data<T_Physic>& update(Type type, const DataStateEvaluator<T_Physic>& ds);
+ };
+};
+
+#endif //H_SPECIFIC_WAVE_DATA_OBJECTIVE_PHASE
diff --git a/specific/model/innerproduct/sobolev.cpp b/specific/model/innerproduct/sobolev.cpp
index 13f09da..7bcfbc9 100644
--- a/specific/model/innerproduct/sobolev.cpp
+++ b/specific/model/innerproduct/sobolev.cpp
@@ -269,11 +269,11 @@ namespace sobolev_df
if(_diag_preconditioner_smooth)
{
- _formulation[c].integral( _weight[c] * _alpha * (_su->get(Order::DIAG,Support::BLK,_mu->get(needToBeUpToDate))[c] + _stabilization[c]) * _D * grad(dof(_j[c])), grad(tf(_j[c])), _config->model_unknown(Support::BLK),integrationType(8*degree));
+ _formulation[c].integral( _weight[c] * _alpha * (_su->get(Order::DIAG,Support::BLK,_mu->get(needToBeUpToDate))[c] + _stabilization[c]) * _D * grad(dof(_j[c])), grad(tf(_j[c])), _config->model_unknown(Support::BLK),integrationType(_integrationDegreeBlk));
}
else
{
- _formulation[c].integral( _weight[c] * _alpha * _p0[c] * _D * grad(dof(_j[c])), grad(tf(_j[c])), _config->model_unknown(Support::BLK),integrationType(8*degree));
+ _formulation[c].integral( _weight[c] * _alpha * _p0[c] * _D * grad(dof(_j[c])), grad(tf(_j[c])), _config->model_unknown(Support::BLK),integrationType(6*degree-2));
}
}
}
@@ -288,7 +288,7 @@ namespace sobolev_df
_filter_formulation.initSystem();
_filter_formulation.integral(dof(_q), tf(_q), _config->model_unknown(Support::BLK), integrationType(2*degree));
- _filter_formulation.integral( _alpha_structure * grad(dof(_q)), grad(tf(_q)), _config->model_unknown(Support::BLK),integrationType(2*degree));
+ _filter_formulation.integral( _alpha_structure * grad(dof(_q)), grad(tf(_q)), _config->model_unknown(Support::BLK),integrationType(2*degree-2));
_filter_formulation.pre();
_filter_formulation.assemble();
@@ -306,7 +306,7 @@ namespace sobolev_df
for (unsigned int c = 0; c < _m_size; c++)
{
//XX
- _filter_formulation.integral(-pow(xComp(grad(m[c])),2), tf(_q), _config->model_unknown(Support::BLK), integrationType(3*degree));
+ _filter_formulation.integral(-pow(xComp(grad(m[c])),2), tf(_q), _config->model_unknown(Support::BLK), integrationType(2*degree-1));
_filter_formulation.assemble();
_filter_formulation.removeTerms();
_filter_formulation.solve(_filterIsFactorized);
@@ -315,7 +315,7 @@ namespace sobolev_df
_m2x[c] = _q;
//YY
- _filter_formulation.integral(-pow(yComp(grad(m[c])),2), tf(_q), _config->model_unknown(Support::BLK), integrationType(3*degree));
+ _filter_formulation.integral(-pow(yComp(grad(m[c])),2), tf(_q), _config->model_unknown(Support::BLK), integrationType(2*degree-1));
_filter_formulation.assemble();
_filter_formulation.removeTerms();
_filter_formulation.solve(_filterIsFactorized);
@@ -324,7 +324,7 @@ namespace sobolev_df
_m2y[c] = _q;
//XY
- _filter_formulation.integral(-xComp(grad(m[c]))*yComp(grad(m[c])), tf(_q), _config->model_unknown(Support::BLK), integrationType(3*degree));
+ _filter_formulation.integral(-xComp(grad(m[c]))*yComp(grad(m[c])), tf(_q), _config->model_unknown(Support::BLK), integrationType(2*degree-1));
_filter_formulation.assemble();
_filter_formulation.removeTerms();
_filter_formulation.solve(_filterIsFactorized);
@@ -380,7 +380,7 @@ namespace sobolev_df
std::complex<double> sum = 0.;
for (unsigned int c = 0; c < _m_size; c++)
{
- sum += _weight[c] * gmshfem::post::integrate( grad(m1[c]) * _D * grad(m2[c]), _config->model_unknown(Support::BLK), integrationType(8*degree) ) / _beta[c] / _beta[c];
+ sum += _weight[c] * gmshfem::post::integrate( grad(m1[c]) * _D * grad(m2[c]), _config->model_unknown(Support::BLK), integrationType(6*degree-2) ) / _beta[c] / _beta[c];
}
return sum;
};
--
GitLab