Skip to content
Snippets Groups Projects
Commit 4db45889 authored by Xavier Adriaens's avatar Xavier Adriaens
Browse files

Phase only objective function

parent 1e1a0e99
No related branches found
No related tags found
No related merge requests found
...@@ -266,6 +266,42 @@ Data<T_Physic> Data<T_Physic>::operator-(const Data<T_Physic>& other) const ...@@ -266,6 +266,42 @@ Data<T_Physic> Data<T_Physic>::operator-(const Data<T_Physic>& other) const
return sub; 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> template<Physic T_Physic>
Data<T_Physic> Data<T_Physic>::operator*(std::complex<double> a) const Data<T_Physic> Data<T_Physic>::operator*(std::complex<double> a) const
{ {
......
...@@ -43,6 +43,8 @@ public: ...@@ -43,6 +43,8 @@ public:
void operator=(const Data<T_Physic>& other); 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*(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;
Data<T_Physic> operator/(std::complex<double> a)const; Data<T_Physic> operator/(std::complex<double> a)const;
private: private:
...@@ -132,18 +134,35 @@ Data<T_Physic> conj(const Data<T_Physic>& d) ...@@ -132,18 +134,35 @@ Data<T_Physic> conj(const Data<T_Physic>& d)
template<Physic T_Physic> template<Physic T_Physic>
Data<T_Physic> abs(const Data<T_Physic>& d) 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 f = 0; f < d.nf(); f++)
{ {
for (unsigned int s = 0; s < d.ns(); s++) for (unsigned int s = 0; s < d.ns(); s++)
{ {
for (unsigned int r = 0; r < d.nr(s); r++) 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 #endif // H_COMMON_DATA_ELEMENT
...@@ -15,6 +15,14 @@ double norm2<Physic::acoustic>(const typename PointData<Physic::acoustic>::Objec ...@@ -15,6 +15,14 @@ double norm2<Physic::acoustic>(const typename PointData<Physic::acoustic>::Objec
template<> template<>
typename PointData<Physic::acoustic>::Object conj<Physic::acoustic>(const typename PointData<Physic::acoustic>::Object& d){return std::conj(d);} typename PointData<Physic::acoustic>::Object conj<Physic::acoustic>(const typename PointData<Physic::acoustic>::Object& d){return std::conj(d);}
template<> 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.;} typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::acoustic>::degree >::Object s_ones<Physic::acoustic>(){return 1.;}
/* electromagnetic */ /* electromagnetic */
...@@ -23,6 +31,17 @@ double norm2<Physic::electromagnetic>(const typename PointData<Physic::electroma ...@@ -23,6 +31,17 @@ double norm2<Physic::electromagnetic>(const typename PointData<Physic::electroma
template<> template<>
typename PointData<Physic::electromagnetic>::Object conj<Physic::electromagnetic>(const typename PointData<Physic::electromagnetic>::Object& d){return d.conjugate();} typename PointData<Physic::electromagnetic>::Object conj<Physic::electromagnetic>(const typename PointData<Physic::electromagnetic>::Object& d){return d.conjugate();}
template<> 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.);} typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::electromagnetic>::degree >::Object s_ones<Physic::electromagnetic>(){return PointData<Physic::electromagnetic>::Object(1.,1.,1.);}
/* elastodynamic */ /* elastodynamic */
...@@ -31,4 +50,15 @@ double norm2<Physic::elastodynamic>(const typename PointData<Physic::elastodynam ...@@ -31,4 +50,15 @@ double norm2<Physic::elastodynamic>(const typename PointData<Physic::elastodynam
template<> template<>
typename PointData<Physic::elastodynamic>::Object conj<Physic::elastodynamic>(const typename PointData<Physic::elastodynamic>::Object& d){return d.conjugate();} typename PointData<Physic::elastodynamic>::Object conj<Physic::elastodynamic>(const typename PointData<Physic::elastodynamic>::Object& d){return d.conjugate();}
template<> 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.);} typename gmshfem::MathObject< std::complex< double >, propertiesOf<Physic::elastodynamic>::degree >::Object s_ones<Physic::elastodynamic>(){return PointData<Physic::elastodynamic>::Object(1.,1.,1.);}
...@@ -28,5 +28,13 @@ template<Physic T_Physic> ...@@ -28,5 +28,13 @@ template<Physic T_Physic>
double norm2(const typename PointData<T_Physic>::Object& d); double norm2(const typename PointData<T_Physic>::Object& d);
template<Physic T_Physic> template<Physic T_Physic>
typename PointData<T_Physic>::Object conj(const typename PointData<T_Physic>::Object& d); 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 #endif // H_SPECIFIC_DATA_ELEMENT
...@@ -9,6 +9,7 @@ ...@@ -9,6 +9,7 @@
//GmshFWI Library //GmshFWI Library
#include "l2distance.h" #include "l2distance.h"
#include "phase.h"
template<Physic T_Physic> 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 = "") 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 ...@@ -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."); throw gmshfem::common::Exception("Objective type could not be found.");
} }
if(objective=="l2distance"){return new l2distance::Objective<T_Physic>(d0,gmshFem);} if(objective=="l2distance"){return new l2distance::Objective<T_Physic>(d0,gmshFem);}
else if(objective=="phase"){return new phase::Objective<T_Physic>(d0,gmshFem);}
else else
{ {
throw gmshfem::common::Exception("Objective " + objective + " is not valid."); throw gmshfem::common::Exception("Objective " + objective + " is not valid.");
......
//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>;
}
#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
...@@ -269,11 +269,11 @@ namespace sobolev_df ...@@ -269,11 +269,11 @@ namespace sobolev_df
if(_diag_preconditioner_smooth) 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 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 ...@@ -288,7 +288,7 @@ namespace sobolev_df
_filter_formulation.initSystem(); _filter_formulation.initSystem();
_filter_formulation.integral(dof(_q), tf(_q), _config->model_unknown(Support::BLK), integrationType(2*degree)); _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.pre();
_filter_formulation.assemble(); _filter_formulation.assemble();
...@@ -306,7 +306,7 @@ namespace sobolev_df ...@@ -306,7 +306,7 @@ namespace sobolev_df
for (unsigned int c = 0; c < _m_size; c++) for (unsigned int c = 0; c < _m_size; c++)
{ {
//XX //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.assemble();
_filter_formulation.removeTerms(); _filter_formulation.removeTerms();
_filter_formulation.solve(_filterIsFactorized); _filter_formulation.solve(_filterIsFactorized);
...@@ -315,7 +315,7 @@ namespace sobolev_df ...@@ -315,7 +315,7 @@ namespace sobolev_df
_m2x[c] = _q; _m2x[c] = _q;
//YY //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.assemble();
_filter_formulation.removeTerms(); _filter_formulation.removeTerms();
_filter_formulation.solve(_filterIsFactorized); _filter_formulation.solve(_filterIsFactorized);
...@@ -324,7 +324,7 @@ namespace sobolev_df ...@@ -324,7 +324,7 @@ namespace sobolev_df
_m2y[c] = _q; _m2y[c] = _q;
//XY //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.assemble();
_filter_formulation.removeTerms(); _filter_formulation.removeTerms();
_filter_formulation.solve(_filterIsFactorized); _filter_formulation.solve(_filterIsFactorized);
...@@ -380,7 +380,7 @@ namespace sobolev_df ...@@ -380,7 +380,7 @@ namespace sobolev_df
std::complex<double> sum = 0.; std::complex<double> sum = 0.;
for (unsigned int c = 0; c < _m_size; c++) 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; return sum;
}; };
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment