diff --git a/common/data/element.cpp b/common/data/element.cpp
index 7195ef09007cfdde4b7349dd2ba38af07ef5111d..3f281fc92c0f2206206804b78f605e992581f5ee 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 cd3667f693900cb4fe737e5dac5592fcb4ce8f18..1198126e131b80d8a0aa29692e9c4611cf83f290 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 53e05138bace5a6a7c37af01e06c37d5dd572e26..03bc76a6187045c2ba7df88ffc1705a00bd69c4c 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 c70de1428c8708da9e173fc5ef56ba877fbf3daf..45654f81c80372b214cbfd47dd0db4bb5b78f5e1 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 0d57e8a4842de7f97cc89b8eaacef623fb89b150..cb6a8a830184de5b039218c3d189a9a03b4800f1 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 0000000000000000000000000000000000000000..8df2b56e6a6a6a80589cb21ba0f0a8f6fd778ea3
--- /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 0000000000000000000000000000000000000000..c653df7ce71630bf2710ee126addfd8c801e0431
--- /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 13f09da4e8a3850af54bf1c9ec04625be954cd4d..7bcfbc9e5de5797417be94178a31a2c8436a2f36 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;
     };