continuous l2 distance up

common/configuration.h

@@ -44,6 +44,8 @@ public:
 
     virtual void mesh() const = 0;
     virtual double area() const = 0;
+    virtual double data_area() const {return gmshfem::post::integrate(1.,_data_omega,"Gauss0");};
+    virtual double datapoint_area() const = 0;
     virtual double array() const {return 0.;};
     virtual double depth() const {return 0.;};
 

common/data/element.cpp

@@ -293,4 +293,32 @@ Data<T_Physic> Data<T_Physic>::phase() const
   return pha;
 }
 
+template<Physic T_Physic>
+Data<T_Physic> Data<T_Physic>::amplitude() const
+{
+  Data<T_Physic> amp(*this);
+  for (unsigned int s = 0; s < ns(); s++)
+  {
+    for (unsigned int r = 0; r < nr(s); r++)
+    {
+      amp.value(s,r, PointData<T_Physic>(std::abs(this->value(s,r)) ) );
+    }
+  }
+  return amp;
+}
+
+template<Physic T_Physic>
+PointData<T_Physic> Data<T_Physic>::minimum() const
+{
+
+  PointData<T_Physic> min(value( (unsigned int) 0,(unsigned int) 0));
+  for (unsigned int s = 0; s < ns(); s++)
+  {
+    for (unsigned int r = 0; r < nr(s); r++)
+    {
+      if(norm(value(s,r)) < norm(min)){min = value(s,r);}
+    }
+  }
+  return min;
+}
 template class Data<Physic::acoustic>;

common/data/element.h

@@ -44,6 +44,8 @@ public:
     Data<T_Physic> operator*(double a)const;
     Data<T_Physic> operator+(double a)const;
     Data<T_Physic> phase()const;
+    Data<T_Physic> amplitude()const;
+    PointData<T_Physic> minimum()const;
 private:
     bool isValid(unsigned int s,unsigned r) const;
 };

input/VO2009Illumination/common_1x1sides.txt

-name = VO2009I_1x1sides
+prename = VO2009I_1x1sides_
 #Configuration
 configuration = inclusion
+receiver_on_emitter = 1
 #
 xe = 0.
 xr = 0.

input/VO2009Illumination/data_1x1sides.txt

 #Data
-data0 = VO2009I_1x1sides_data0
-data1 = VO2009I_1x1sides_data1
-data2 = VO2009I_1x1sides_data2
-data3 = VO2009I_1x1sides_data3
-data4 = VO2009I_1x1sides_data4
-data5 = VO2009I_1x1sides_data5
+data0 = VO2009I_1x1sides_synthetics_data0
+data1 = VO2009I_1x1sides_synthetics_data1
+data2 = VO2009I_1x1sides_synthetics_data2
+data3 = VO2009I_1x1sides_synthetics_data3
+data4 = VO2009I_1x1sides_synthetics_data4
+data5 = VO2009I_1x1sides_synthetics_data5

input/VO2009Illumination/gradient.txt

+name = gradient
+#Configuration
+inclusion = none
+unknown = background
+#
+#Objective
+objective = l2distance
+#
+#InnerProduct
+diag_preconditioner = 0
+ref_preconditioner = 1.3e5
+complex = 0
+#Scale
+scale0 = 10.
+scaleN = 10.
+interval = linear
+N = 0
+#
+#Discretization
+h = 0.05
+model_FunctionSpaceType = HierarchicalH1
+model_FunctionSpaceDegree = 5

input/VO2009Illumination/inversion.txt

+name = inversion
 #Configuration
 inclusion = none
 unknown = background
@@ -8,9 +9,11 @@ gaussNewton = 1
 objective = l2distance
 #
 #Inner product
-complex = 0
 innerproduct = sobolev
-diag_preconditioner = 1
+boundaryWeight = 0.
+diag_preconditioner = 0
+ref_preconditioner = 1e3
+complex = 0
 #
 n_group=6
 #Discretization
@@ -20,14 +23,32 @@ model_IntegrationType = Gauss
 model_FunctionSpaceType = HierarchicalH1
 model_FunctionSpaceDegree = 1
 #
-#Minimum search
-minimumsearch = classic
-minimum_maxIteration = 2
+#Global minimum search
+globalminimumsearch = sequentialscale
+globalminimum_n_scale = 1
+globalminimum_minNscale = 1
+globalminimum_relDecreaseThreshold = 0.
+#
+scale00 = 1.
+scale10 = 1.
+scale20 = 1.
+scale30 = 1.
+scale40 = 1.
+scale50 = 1.
+#
+#Local minimum search
+localminimumsearch = classic
+#localminimum_writeInterval = 5
+localminimum_maxIteration = 4
+#localminimum_absDecreaseThreshold = 0.
+#localminimum_overallRelDecreaseThreshold = 1.
+#localminimum_meanRelDecreaseThreshold = 0.
+#localminimum_meanRelDecreaseMemory = 22
 #
 #Descent search
 descentsearch = gradientdescent
 #
 #Line search
-linesearch = backtracking
+linesearch = armijo
 rho = 0.5
-c = 0.3
+c1 = 0.3

input/VO2009Illumination/preconditioner.txt

+name = preconditioner
+#Configuration
+inclusion = none
+unknown = background
+#
+#Objective
+objective = l2distance
+#Discretization
+h = 0.05
+#Model
+model_FunctionSpaceDegree = 1

input/VO2009Illumination/synthetics.txt

+name = synthetics
 #Configuration
 inclusion = cylinder
 filled = 1
@@ -10,5 +11,7 @@ Im(mc) = 0.
 #
 #Discretization
 h = 0.05
+model_FunctionSpaceDegree = 1
 #Output
-write_fields = 0
+write_wave_fields = 0
+write_data_fields = 0

input/paper3/common.txt

+prename = paper2_
+#Configuration
+configuration = soil
+#
+xoffset = 0
+Ler = 8.712
+nr = 122
+ne = 122
+ReceiverOnEmitter=1
+ymax = 0.216
+ymin = 0.072
+L = 9.192
+H = 2.904
+#
+Re(v_air) = 1.5
+Im(v_air) = 0.
+Re(v_soil) = 2.
+Im(v_soil) = 0.
+#
+#Equation
+physic = acoustic
+equation = helmholtz
+#
+#Discretization
+h = 0.036
+#Wave
+wave_IntegrationType = Gauss
+wave_FunctionSpaceType = HierarchicalH1
+#
+#Frequencies
+n_freq = 3
+#Frequency0
+frequency0 = 4.
+#Frequency1
+frequency1 = 6.
+#Frequency2
+frequency2 = 8.

input/paper3/inversion.txt

+#Configuration
+unknown = soil
+m0_type = velocity_depth
+Re(v_0) = 1.5
+Im(v_0) = 0.
+Re(v_H) = 3.5
+Im(v_H) = 0.
+#
+#Equation
+gaussNewton=0
+#
+#Global minimum search
+globalminimum_n_scale=1
+#
+scale00 = 1.5
+#
+scale10 = 0.35
+#
+scale20  = 0.25
+#
+n_group=3
+#Data
+data0 = paper2_synthetics_data0
+data1 = paper2_synthetics_data1
+data2 = paper2_synthetics_data2
+#
+#Discretization
+#Wave
+wave_FunctionSpaceDegree0 = 2
+wave_FunctionSpaceDegree1 = 3
+wave_FunctionSpaceDegree2 = 4
+#Model
+model_IntegrationType = Gauss
+model_FunctionSpaceType = HierarchicalH1
+model_FunctionSpaceDegree0 = 1
+model_FunctionSpaceDegree1 = 1
+model_FunctionSpaceDegree2 = 1
+#
+#Global minimum search
+globalminimumsearch = sequentialscale
+globalminimum_relDecreaseThreshold = 0.0
+#Local minimum search
+localminimum_absDecreaseThreshold = 1e-5
+#
+#Inner product
+complex = 0
+boundaryWeight = 0.
+diag_preconditioner = 1
+innerproduct = sobolev

input/paper3/inversion_continuous.txt

+name = inversion_continuous
+#
+#Objective
+objective = continuousl2distance
+objective_l_s = 1.5
+objective_weight = paper2_reference_data0
+#Inner product
+ref_preconditioner = 50
+#Local minimum search
+localminimumsearch = yuanfan
+localminimum_initMu = 1
+localminimum_c0 = 0.0001
+localminimum_c1 = 0.75
+localminimum_c2 = 0.75
+localminimum_c3 = 0.25
+localminimum_c4 = 2.
+localminimum_writeInterval = 5
+localminimum_maxIteration = 30
+#
+#Descent search
+descentsearch = newton_steihaug_conjugategradient
+descent_maxIteration = 50
+descent_eta0 = 0.4

input/paper3/ip_comparison.txt

+name = ip_comparison
+#Configuration
+unknown = soil
+m0_type = velocity_depth
+Re(v_0) = 1.5
+Im(v_0) = 0.
+Re(v_H) = 3.5
+Im(v_H) = 0.
+#
+#Equation
+gaussNewton=0
+#Objective
+objective = continuousl2distance
+objective_weight = paper2_reference_data0
+objective_l_s = 1.5
+#
+#Inner product
+complex = 0
+boundaryWeight = 0.
+diag_preconditioner = 1
+innerproduct = sobolev
+ref_preconditioner = 50
+#
+scale0 = 1.5
+#
+#Data
+data0 = paper2_synthetics_data0
+#
+#Discretization
+#Wave
+wave_FunctionSpaceDegree0 = 2
+#Model
+model_IntegrationType = Gauss
+model_FunctionSpaceType = HierarchicalH1
+model_FunctionSpaceDegree = 1

input/paper3/reference.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Thu May  6 15:15:24 2021
+
+@author: xavier
+"""
+import numpy as np
+import scipy as sp
+import scipy.special
+import matplotlib.pyplot as plt
+from matplotlib import rc
+rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb}')
+rc('font',**{'family':'serif','sans-serif':['Helvetica'],'size':10})
+rc('text', usetex=True)
+#Value for 1 width figure
+H = 1.5;
+W = 6.5;
+save = False;
+
+
+"""
+Example I: Wave number perturbation
+"""
+L = 8.712
+f = 4.
+ω = 2. * np.pi * f
+v = 1.5
+k = ω / v
+N = 122
+x = np.linspace(0, L, N)
+
+def d(xe,xr):
+        if(xe==xr): return 0.5 * (d(0.,1e-3)+d(0.,-1e-3)) #Must match numerical value at 0 (?)
+        return np.conj( -1j/4*sp.special.hankel1(0,k*np.abs(xe-xr)) )
+
+D = np.zeros([N, 2*N])
+for e in range(0,N):
+    for r in range(0,N):
+        D[e,2*r] = np.real(d(x[e],x[r]))
+        D[e,2*r+1] = np.imag(d(x[e],x[r]))
+
+np.savetxt("paper2_reference_data0.csv",D,delimiter=";")
\ No newline at end of file

input/paper3/synthetics.txt

+name = synthetics
+#Configuration
+unknown = none
+m0_type = velocity_file
+path = ../input/paper2/marmousi
+#
+#Discretization
+#Wave
+wave_FunctionSpaceDegree0 = 2
+wave_FunctionSpaceDegree1 = 3
+wave_FunctionSpaceDegree2 = 4
+#Model
+model_IntegrationType = Gauss
+model_FunctionSpaceDegree = 4
+#Output
+write_data_fields = 1
+write_wave_fields = 0
+#write_single_wave_fields = 3
+#write_single_wave_field0 = 30
+#write_single_wave_field1 = 61
+#write_single_wave_field2 = 92

specific/configuration/inclusion.cpp

@@ -620,10 +620,10 @@ namespace inclusion
             int cl2 = factory::addCurveLoop({lyHr2,lxLe,lyLr,lxHe2});
             int s2 = factory::addPlaneSurface({cl2});
 
-            factory::mesh::setTransfiniteCurve(lyHr2,_nxe);
-            factory::mesh::setTransfiniteCurve(lxLe,_nxr);
-            factory::mesh::setTransfiniteCurve(lyLr,_nxe);
-            factory::mesh::setTransfiniteCurve(lxHe2,_nxr);
+            factory::mesh::setTransfiniteCurve(lyHr2,_nye);
+            factory::mesh::setTransfiniteCurve(lxLe,_nyr);
+            factory::mesh::setTransfiniteCurve(lyLr,_nye);
+            factory::mesh::setTransfiniteCurve(lxHe2,_nyr);
 
             factory::mesh::setTransfiniteSurface(s2,"Left",{pHeHr,pLe0,pLeLr,p0Lr});
 
@@ -700,16 +700,22 @@ namespace inclusion
     }
 
     bool Configuration::data_coordinate_isValid(double xs,double xr) const
+    {
+        double eps = 1e-8;
+        gmsh::option::getNumber("Geometry.Tolerance",eps);
+        if( nxe() != 0 )
         {
             double xe0 = (_H-_He)/2.;
             double xr0 = (_H-_Hr)/2.;
-        if(xe0 <= xs && xs <= xe0+_He && xr0 <= xr && xr <= xr0+_Hr){return true;}
-
-        xe0 = _H + (_L-_Le)/2.;
-        xr0 = _H + (_L-_Lr)/2.;
-        if(xe0 <= xs && xs <= xe0+_He && xr0 <= xr && xr <= xr0+_Hr){return true;}
+            if( (xe0-eps <= xs) && (xs <= xe0+_He+eps) && (xr0-eps <= xr) && (xr <= xr0+_Hr+eps) ){return true;}
+        }
 
-        throw Exception("data_coordinate_isValid not implemented for inclusion");
+        if( nye() != 0 )
+        {
+            double xe0 = _H + (_L-_Le)/2.;
+            double xr0 = _H + (_L-_Lr)/2.;
+            if( (xe0-eps <= xs) && (xs <= xe0+_Le+eps) && (xr0-eps <= xr) && (xr <= xr0+_Lr+eps) ){return true;}
+        }
         return false;
     }
 

specific/configuration/inclusion.h

@@ -102,6 +102,7 @@ namespace inclusion
         virtual bool data_coordinate_isValid(double xs,double xr) const;
 
         virtual double area() const {return _L * _H; };
+        virtual double datapoint_area() const {return _Le * _Lr / (double(_nye)) / (double(_nyr)) + _He * _Hr / (double(_nxe)) / (double(_nxr));};
         virtual double array() const {return std::max(_He,_Le);};
         virtual double depth() const {return std::max(_H,_L)/2.;};
         virtual double m_reference() const {return std::real(_mb);};

specific/configuration/soil.cpp

@@ -12,6 +12,7 @@
 #include "gmsh.h"
 //Fwi library
 #include "../../common/wave/correlation.h"
+#include "../../common/data/element.h"
 #include "soil.h"
 
 namespace gmodel = gmsh::model;
@@ -519,20 +520,28 @@ namespace soil
 
     bool Configuration::data_coordinate_isValid(double xs,double xr) const
     {
-        //Bof bof correct !!
-        bool valid = (0. <= xs) && (xs <= _L) && (0. <= xr) && (xr <= _L);
-        if(!_receiverOnEmitter)
-        {
-            valid = valid && !(std::abs(xr-xs) < 1e-6 * std::abs(xr+xs));
-        }
-        return valid;
+        double eps = 1e-8;
+        gmsh::option::getNumber("Geometry.Tolerance",eps);
+        if( (!_receiverOnEmitter) && (std::abs(xr-xs) < eps) ){return false;}
+
+        double x0 = (_L-_Ler)/2.;
+        return (x0-eps <= xs) && (xs <= x0+_Ler+eps) && (x0-eps <= xr) && (xr <= x0+_Ler+eps);
     }
 
     template<Physic T_Physic>
     ModelFunction green0_preconditioner(const Data<T_Physic>& dd, const WaveMultiField<T_Physic>& g, const Configuration* const config)
     {
+        if(dd.constant())
+        {
+            unsigned int e = 0;
+            unsigned int r = 0;
+            return dd.value(e,r) * pow(autocorrelate<T_Physic>(g,config->all_emitter_idx(true)),2);
+        }
+        else
+        {
+            return pow(autocorrelate<T_Physic>(g,config->all_emitter_idx(true)),2);
             throw common::Exception("green0_preconditioner not implemented for soil");
-        //return pow(autocorrelate<T_Physic>(g,config->all_emitter_idx(true)),2);
+        }
     }
 
     template ModelFunction green0_preconditioner<Physic::acoustic>(const Data<Physic::acoustic>&, const WaveMultiField<Physic::acoustic>&, const Configuration* const);

specific/configuration/soil.h

@@ -60,6 +60,7 @@ namespace soil
 
     double H() const {return _H;};
     virtual double area() const {return _L * _H;};
+    virtual double datapoint_area() const {return _Ler * _Ler / (double(_ne)) / (double(_nr));};
     virtual double array() const {return _Ler;};
     virtual double depth() const {return _H;};
     virtual double m_reference() const {return std::real(_m_soil);};

specific/data/objective/continuousl2distance.cpp

+//GmshFEM Library
+#include "Exception.h"
+#include "Function.h"
+#include "AlgebraicFunctions.h"
+
+//GmshFWI Library
+#include "continuousl2distance.h"
+
+#define d0 *_d0
+#define d ds.state(Type::F)
+#define rho ds.state(Type::A)
+#define dd ds.state(Type::PF)
+#define drho ds.state(Type::PA)
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+using namespace gmshfem::function;
+using namespace gmshfem::equation;
+using namespace gmshfem::post;
+/*
+* SOBOLEV
+*/
+namespace continuousl2distance
+{
+    template<Physic T_Physic>
+    void  Objective<T_Physic>::updateSystem_s()
+    {
+        if(_systemIsUpToDate_s){return;}
+        _formulation_s.removeSystem();
+        _formulation_s.initSystem();
+
+        _formulation_s.galerkin( dof(_v_s), tf(_v_s), _config->data_omega(), "Gauss2");
+
+        TensorFunction< std::complex< double > > Cs = (_l_s) / 2. / M_PI * tensor< std::complex< double > >(1., 0., 0., 0., 1., 0., 0., 0., 0.);
+        _formulation_s.galerkin(  Cs * grad(dof(_v_s)), Cs * grad(tf(_v_s)), _config->data_omega(), "Gauss2");
+
+        _formulation_s.pre();
+        _formulation_s.assemble();
+        _systemIsUpToDate_s=true;
+        _formulation_s.removeTerms();
+    }
+
+    template<Physic T_Physic>
+    Data<T_Physic> Objective<T_Physic>::inv_operator_s(const Data<T_Physic>& r)
+    {
+        DataField<T_Physic> rF(_v_s);
+        rF.setValuesToZero();
+        rF.forward_interpolate(r);
+
+        if(!_systemIsUpToDate_s){updateSystem_s();}
+        _formulation_s.galerkin( -_a_s * rF, tf(_v_s), _config->data_omega(), "Gauss2");
+        _formulation_s.assemble();
+        _formulation_s.removeTerms();
+
+        _formulation_s.solve(_systemIsFactorized_s);
+        _systemIsFactorized_s=true;
+        _formulation_s.setRHSToZero();
+        rF = _v_s;
+
+        Data<T_Physic> out(r);
+        out.value(rF);
+        return out;
+    }
+
+    template<Physic T_Physic>
+    double Objective<T_Physic>::performance(const Data<T_Physic>& d1)
+    {
+
+        Data<T_Physic> r = (d1 - d0) / _w;
+        ( d1 ).write("d1","pos");
+        ( d0 ).write("d0","pos");
+        ( d1 / _w ).write("d1w","pos");
+        ( d0 / _w ).write("d0w","pos");
+        r.write("r","pos");
+        Data<T_Physic> sd1w = inv_operator_s( d1 / _w );
+        Data<T_Physic> sd0w = inv_operator_s( d0 / _w );
+        Data<T_Physic> sr = inv_operator_s(r);
+        sd0w.write("sd0w","pos");
+        sd1w.write("sd1w","pos");
+        (sd1w.phase() / sd0w.phase()).write("sd01w","pos");
+        sr.write("sr","pos");
+        double coef =  _config->datapoint_area() / _config->data_area();
+        //msg::print << "coef = " << coef << msg::endl;
+        return 0.5 * coef * std::real( l2scalp(sr,sr) );
+        //return 0.5 * std::real( integrate(_v_s * conj(_v_s), _config->data_omega(), "Gauss2") );
+    }
+
+    template<Physic T_Physic>
+    const Data<T_Physic>& Objective<T_Physic>::update(Type type, const DataStateEvaluator<T_Physic>& ds)
+    {
+        double coef =  _config->datapoint_area() / _config->data_area();
+        //msg::print << "coef = " << coef << msg::endl;
+        switch (type)
+        {
+            case Type::F: case Type::PF: default:
+                throw Exception("Objective can not update (perturbed) forward data.");
+                break;
+            case Type::A:
+                _v = conj( inv_operator_s((inv_operator_s((d-d0)/_w))) ) / _w * coef;
+                _v.write("_v","pos");
+                break;
+            case Type::PA:
+                _v = conj( inv_operator_s((inv_operator_s((dd)/_w))) ) / _w * coef;
+                break;
+            case Type::D:
+                _v.value(PointData<T_Physic>(coef / norm(_w.minimum()) / norm(_w.minimum()) ));
+                break;
+        }
+        return _v;
+    }
+    template class Objective<Physic::acoustic>;
+}