reproducing Figure 1 of Virieux & Operto 2009 review paper

CONVENTIONS.txt

 /*
 * Content
 */
-*Adjoint/Dual fields are always stored through their conjugate
+*Adjoint/Dual fields are always stored through their conjugate for data, wave and sensitivities but not for model.
 
 
 /*

common/configuration.h

@@ -21,6 +21,8 @@
 class ConfigurationInterface //(abstract)
 {
 protected:
+    const std::string _name;
+
     std::array<gmshfem::domain::Domain,2> _wave_omega;
     std::array<gmshfem::domain::Domain,2> _model_known;
     std::array<gmshfem::domain::Domain,2> _model_unknown;
@@ -40,7 +42,7 @@ protected:
     virtual void wave_mesh() const = 0;
     virtual void data_mesh() const = 0;
 public:
-    ConfigurationInterface(const gmshfem::common::GmshFem& gmshFem)
+    ConfigurationInterface(std::string name, const gmshfem::common::GmshFem& gmshFem) : _name(name)
     {
         unsigned int remesh = true;
         gmshFem.userDefinedParameter(remesh, "remesh");
@@ -49,6 +51,7 @@ public:
     virtual ~ConfigurationInterface() = default;
 
     virtual void mesh() const = 0;
+    virtual double area() const = 0;
     virtual double array() const {return 0.;};
     virtual double depth() const {return 0.;};
 

common/functional.h

@@ -3,12 +3,15 @@
 
 //GmshFWI Library
 #include "model/element.h"
+#include "model/innerproduct/innerproduct.h"
 
 class FunctionalInterface
 {
 public:
     virtual ~FunctionalInterface() = default;
 
+    virtual const InnerProductInterface* const innerproduct() const = 0;
+
     virtual void setModel(const ModelField& m) = 0;
     virtual void setModel(const ModelFunction& m) = 0;
     virtual void setModelPerturbation(const ModelField& dm) = 0;

common/model/discretization.cpp

@@ -14,14 +14,23 @@ namespace model
   /*
   *  class Discretization
   */
-  Discretization::Discretization(const GmshFem& gmshFem)
+  Discretization::Discretization(const GmshFem& gmshFem, std::string suffix)
   {
     msg::print << "Read model discretization parameters" << msg::endl;
     std::string str_buffer;
-    if(
-      !(gmshFem.userDefinedParameter(_functionSpaceDegree, "model_FunctionSpaceDegree") &&
-      gmshFem.userDefinedParameter(str_buffer, "model_FunctionSpaceType"))
+    if(!(
+            (
+                gmshFem.userDefinedParameter(_functionSpaceDegree, "model_FunctionSpaceDegree")
+                ||
+                gmshFem.userDefinedParameter(_functionSpaceDegree, "model_FunctionSpaceDegree"+suffix)
             )
+            &&
+            (
+                gmshFem.userDefinedParameter(str_buffer, "model_FunctionSpaceType")
+                ||
+                gmshFem.userDefinedParameter(str_buffer, "model_FunctionSpaceType"+suffix)
+            )
+    ))
     {
       throw common::Exception("A model discretization parameter could not be found.");
     }

common/model/discretization.h

@@ -16,7 +16,7 @@ namespace model
       gmshfem::field::FunctionSpaceTypeForm0 _functionSpaceType;
       unsigned int _functionSpaceDegree;
   public:
-      Discretization(const gmshfem::common::GmshFem& gmshFem);
+      Discretization(const gmshfem::common::GmshFem& gmshFem, std::string suffix = "");
       gmshfem::field::FunctionSpaceTypeForm0 functionSpaceType() const {return _functionSpaceType;};
       unsigned int functionSpaceDegree() const {return _functionSpaceDegree;};
       void functionSpaceDegree(unsigned int order) {_functionSpaceDegree=order;};

common/model/element.cpp

@@ -11,6 +11,25 @@ using namespace gmshfem;
 using namespace gmshfem::common;
 using namespace gmshfem::post;
 
+/*
+* ModelField
+*/
+void ModelField::oppose()
+{
+    for (auto it = this->begin(); it != this->end(); it++)
+    {
+        it->second = -it->second;
+    }
+}
+
+void ModelField::conjugate()
+{
+    for (auto it = this->begin(); it != this->end(); it++)
+    {
+        it->second = std::conj(it->second);
+    }
+}
+
 /*
 * Model
 */

common/model/element.h

@@ -20,6 +20,9 @@ public:
   ModelField(const gmshfem::field::Field< std::complex<double> , gmshfem::field::Form::Form0 >& other) : gmshfem::field::Field< std::complex<double> , gmshfem::field::Form::Form0 >(other) {};
 
   void write(std::string name, std::string type="pos", std::string path="") const {gmshfem::post::save(*this,this->domain(),name,type,path);}
+
+  void oppose();
+  void conjugate();
 };
 
 /*

common/model/innerproduct/innerproduct.cpp

@@ -13,10 +13,14 @@ static const std::complex< double > im = std::complex< double >(0., 1.);
 /*
 * DifferentialInnerProductInterface
 */
-DifferentialInnerProductInterface::DifferentialInnerProductInterface(bool complex, const ConfigurationInterface* const config, const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem)
+DifferentialInnerProductInterface::DifferentialInnerProductInterface(bool complex, const ConfigurationInterface* const config, const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem,std::string suffix)
 : InnerProductInterface(complex, config, m_discret), _j("_j",_config->model_evaldomain(),_config->model_gmodel(),_m_discret), _formulation("model_formulation"),_systemIsUpToDate(false),_systemIsFactorized(false)
 {
-    if(!gmshFem.userDefinedParameter(_integrationType, "model_IntegrationType"))
+    if(!(
+        gmshFem.userDefinedParameter(_integrationType, "model_IntegrationType")
+        ||
+        gmshFem.userDefinedParameter(_integrationType, "model_IntegrationType"+suffix)
+    ))
     {
       throw common::Exception("Model integration type could not be found.");
     }
@@ -56,7 +60,7 @@ const ModelField& DifferentialInnerProductInterface::update(const ModelFunction&
     return _j;
 }
 
-std::complex<double> DifferentialInnerProductInterface::product(const ModelField& m1, const ModelField& m2)
+std::complex<double> DifferentialInnerProductInterface::product(const ModelField& m1, const ModelField& m2) const
 {
     algebra::MatrixCRSFast< std::complex<double> > M;
     algebra::Vector< std::complex<double> > m1v, m2v;

common/model/innerproduct/innerproduct.h

@@ -31,7 +31,7 @@ public:
     virtual void unlink() {};
     virtual const ModelField& update(const ModelFunction& bulk_sensi, const ModelFunction& boundary_sensi) = 0;
 
-    virtual std::complex<double> product(const ModelField& m1, const ModelField& m2) = 0;
+    virtual std::complex<double> product(const ModelField& m1, const ModelField& m2) const = 0;
 
     virtual void modelIsObsolete() {};
 };
@@ -49,11 +49,11 @@ protected:
     bool _systemIsUpToDate;
     bool _systemIsFactorized;
 public:
-    DifferentialInnerProductInterface(bool complex, const ConfigurationInterface* const config, const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem);
+    DifferentialInnerProductInterface(bool complex, const ConfigurationInterface* const config, const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem, std::string suffix = "" );
     virtual const ModelField& update(const ModelFunction& bulk_sensi, const ModelFunction& boundary_sensi) ;
     std::string integrationType(unsigned int degree) const {return _integrationType + std::to_string(degree);};
 
-    virtual std::complex<double> product(const ModelField& m1, const ModelField& m2);
+    virtual std::complex<double> product(const ModelField& m1, const ModelField& m2) const;
 protected:
     virtual void setLHS() = 0;
     virtual void setRHS(const ModelFunction& bulk_sensi, const ModelFunction& boundary_sensi) = 0;

common/optimization/descentsearch.h

@@ -3,6 +3,7 @@
 
 //GmshFWI Library
 #include "../functional.h"
+#include "../model/innerproduct/innerproduct.h"
 
 /*
 * class DescentSearchHistoryInterface
@@ -25,7 +26,7 @@ class DescentSearchInterface
 public:
     virtual ~DescentSearchInterface() = default;
 
-    virtual void descentsearch(FunctionalInterface* const functional) const = 0;
+    virtual void operator()(FunctionalInterface* const functional) const = 0;
 
     virtual const DescentSearchHistoryInterface* const history() const = 0;
 };

common/optimization/linesearch.h

@@ -28,7 +28,7 @@ class LineSearchInterface
 public:
     virtual ~LineSearchInterface() = default;
 
-    virtual double linesearch(FunctionalInterface* const functional) const = 0;
+    virtual double operator()(FunctionalInterface* const functional) const = 0;
 
     virtual const LineSearchHistoryInterface* const history() const = 0;
 };

common/optimization/minimumsearch.h

@@ -24,7 +24,7 @@ class MinimumSearchInterface
 public:
     virtual ~MinimumSearchInterface() = default;
 
-    virtual void minimize(ModelField* const m, FunctionalInterface* const functional, const DescentSearchInterface* const descentsearch, const LineSearchInterface* const linesearch) const = 0;
+    virtual void operator()(ModelField* const m, FunctionalInterface* const functional, const DescentSearchInterface* const descentsearch, const LineSearchInterface* const linesearch) const = 0;
 
     virtual const MinimumSearchHistoryInterface* const history() const = 0;
 };

common/rediscretize.cpp

+//GmshFEM Library
+#include "Formulation.h"
+#include "Function.h"
+
+//GmshFWI Library
+#include "rediscretize.h"
+
+using namespace gmshfem::problem;
+using namespace gmshfem::equation;
+using namespace gmshfem::function;
+
+ModelField rediscretize(const ModelField& m0, const ConfigurationInterface* const config, const model::Discretization& m_discret, std::string integrationType)
+{
+    integrationType += std::to_string(2*m_discret.functionSpaceDegree());
+    ModelField m("m",config->model_evaldomain(), config->model_gmodel(),m_discret);
+    Formulation< std::complex< double > > formulation("Model rediscretization");
+    formulation.galerkin(dof(m), tf(m),config->model_unknown(Support::BLK), integrationType);
+    formulation.galerkin(-1. * m0, tf(m),config->model_unknown(Support::BLK),integrationType);
+    formulation.pre();
+    formulation.assemble();
+    formulation.solve();
+    return m;
+}

common/rediscretize.h

+#ifndef H_COMMON_REDISCRETIZE
+#define H_COMMON_REDISCRETIZE
+
+//GmshFWI
+#include "model/element.h"
+#include "configuration.h"
+
+ModelField rediscretize(const ModelField & m0, const ConfigurationInterface* const config, const model::Discretization& m_discret, std::string integrationType);
+
+#endif // H_COMMON_REDISCRETIZE

common/statefunctional.h

@@ -36,6 +36,8 @@ public:
     void write_data(Type type, std::string name);
     void write_wave(Type type, std::string name);
 
+    virtual const InnerProductInterface* const innerproduct() const {return _innerproduct;};
+
     virtual void setModel(const ModelField& m);
     virtual void setModel(const ModelFunction& m);
     virtual void setModelPerturbation(const ModelField& dm);

convergence.cpp

@@ -35,12 +35,12 @@ int convergence(const GmshFem& gmshFem)
     std::string name = "noname";
     gmshFem.userDefinedParameter(name, "name");
 
-    ConfigurationInterface* configuration = newConfiguration(gmshFem);
+    ConfigurationInterface* configuration = newConfiguration(name, gmshFem);
 
     wave::Discretization<T_Physic> w_discret(gmshFem);
 
     double frequency = 0.;
-    if(!gmshFem.userDefinedParameter(frequency, "frequency"))
+    if(!gmshFem.userDefinedParameter(frequency, "frequency0"))
     {
         throw Exception("Frequency could not be found.");
     }
@@ -48,7 +48,7 @@ int convergence(const GmshFem& gmshFem)
     EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem);
 
     std::string filename = "noname";
-    gmshFem.userDefinedParameter(filename, "data");
+    gmshFem.userDefinedParameter(filename, "data0");
     Data<T_Physic> d0(filename,configuration);
 
     ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(&d0,gmshFem);

directional.cpp

@@ -35,12 +35,12 @@ int directional(const GmshFem& gmshFem)
     std::string name = "noname";
     gmshFem.userDefinedParameter(name, "name");
 
-    ConfigurationInterface* configuration = newConfiguration(gmshFem);
+    ConfigurationInterface* configuration = newConfiguration(name, gmshFem);
 
     wave::Discretization<T_Physic> w_discret(gmshFem);
 
     double frequency = 0.;
-    if(!gmshFem.userDefinedParameter(frequency, "frequency"))
+    if(!gmshFem.userDefinedParameter(frequency, "frequency0"))
     {
         throw Exception("Frequency could not be found.");
     }
@@ -48,7 +48,7 @@ int directional(const GmshFem& gmshFem)
     EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem);
 
     std::string filename = "noname";
-    gmshFem.userDefinedParameter(filename, "data");
+    gmshFem.userDefinedParameter(filename, "data0");
     Data<T_Physic> d0(filename,configuration);
 
     ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(&d0,gmshFem);

gradient.cpp

@@ -36,11 +36,11 @@ int gradient(const GmshFem& gmshFem)
     std::string name = "noname";
     gmshFem.userDefinedParameter(name, "name");
 
-    ConfigurationInterface* configuration = newConfiguration(gmshFem);
+    ConfigurationInterface* configuration = newConfiguration(name, gmshFem);
 
 
     double frequency = 0.;
-    if(!gmshFem.userDefinedParameter(frequency, "frequency"))
+    if(!gmshFem.userDefinedParameter(frequency, "frequency0"))
     {
         throw Exception("Frequency could not be found.");
     }
@@ -49,7 +49,7 @@ int gradient(const GmshFem& gmshFem)
     EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem);
 
     std::string filename = "noname";
-    gmshFem.userDefinedParameter(filename, "data");
+    gmshFem.userDefinedParameter(filename, "data0");
     Data<T_Physic> d0(filename,configuration);
 
     ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(&d0,gmshFem);

input/VO2009Algorithm/common.txt

+#Configuration
+configuration = inclusion
+#
+xe = 0.2
+xr = 2.8
+ye = 0.2
+yr = 2.8
+He = 2.6
+Hr = 2.6
+Le = 2.6
+Lr = 2.6
+nxe = 14
+nxr = 14
+nye = 14
+nyr = 14
+L = 3.0
+H = 3.0
+#
+Re(mb) = 1.
+Im(mb) = 0.
+#
+#Equation
+physic = acoustic
+equation = helmholtz
+#
+#Discretization
+#Wave
+wave_IntegrationType = Gauss
+wave_FunctionSpaceType = HierarchicalH1
+#
+n_freq = 4
+#Frequency0
+frequency0 = 1.
+wave_FunctionSpaceDegree0 = 2
+#
+#Frequency1
+frequency1 = 1.25
+wave_FunctionSpaceDegree1 = 2
+#
+#Frequency2
+frequency2 = 1.75
+wave_FunctionSpaceDegree2 = 3
+#
+#Frequency3
+frequency3 = 2.5
+wave_FunctionSpaceDegree3 = 3

input/VO2009Algorithm/inversion_gd.txt

+name = VO2009A_gd
+#Configuration
+inclusion = none
+unknown = background
+#
+#Data
+data0 = VO2009A_data0
+data1 = VO2009A_data1
+data2 = VO2009A_data2
+data3 = VO2009A_data3
+#
+#Equation
+gaussNewton = 1
+#Objective
+objective = l2distance
+#
+#Inner product
+complex = 0
+innerproduct = sobolev
+diag_preconditioner = 1
+#
+n_group=4
+#Discretization
+h = 0.1
+#Model
+model_IntegrationType = Gauss
+model_FunctionSpaceType = HierarchicalH1
+model_FunctionSpaceDegree = 1
+#
+#Minimum search
+minimumsearch = classic
+minimum_maxIteration = 1
+#
+#Descent search
+descentsearch = gradientdescent
+descent_maxIteration = 1
+#
+#Line search
+linesearch = secondorder