Minimum search

CMakeLists.txt

@@ -111,3 +111,5 @@ add_executable( gradient gradient.cpp ${EXTRA_INCS})
 target_link_libraries( gradient fwi  ${EXTRA_LIBS})
 add_executable( convergence convergence.cpp ${EXTRA_INCS})
 target_link_libraries( convergence fwi  ${EXTRA_LIBS})
+add_executable( inversion inversion.cpp ${EXTRA_INCS})
+target_link_libraries( inversion fwi  ${EXTRA_LIBS})

CONVENTIONS.txt

+/*
+* Content
+*/
+*Adjoint/Dual fields are always stored through their conjugate
+
+
+/*
+* Implementation
+*/
+*Passing through references is used only to avoid copy. Passing through pointers is used otherwise.  

common/initialization.cpp

+//GmshFEM Library
+#include "Formulation.h"
+
+//GmshFWI Library
+#include "initialization.h"
+
+using namespace gmshfem::problem;
+using namespace gmshfem::equation;
+
+ModelField initialize(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 initialization");
+    formulation.galerkin(dof(m), tf(m),config->model_unknown(Support::BLK), integrationType);
+    formulation.galerkin(-config->m0(), tf(m),config->model_unknown(Support::BLK),integrationType);
+    formulation.pre();
+    formulation.assemble();
+    formulation.solve();
+    return m;
+}

common/initialization.h

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

common/model/discretization.h

@@ -15,15 +15,9 @@ namespace model
   private:
       gmshfem::field::FunctionSpaceTypeForm0 _functionSpaceType;
       unsigned int _functionSpaceDegree;
-      std::string _integrationType;
   public:
       Discretization(const gmshfem::common::GmshFem& gmshFem);
       gmshfem::field::FunctionSpaceTypeForm0 functionSpaceType() const {return _functionSpaceType;};
-      std::string integrationType(unsigned int integrandFieldNumber) const
-      {
-          return _integrationType + std::to_string(integrandFieldNumber*_functionSpaceDegree);
-      };
-      std::string integrationType() const {return _integrationType;};
       unsigned int functionSpaceDegree() const {return _functionSpaceDegree;};
       void functionSpaceDegree(unsigned int order) {_functionSpaceDegree=order;};
   };

common/model/element.h

@@ -18,6 +18,7 @@ public:
   ModelField(std::string name,const gmshfem::domain::Domain &domain, std::string gmodel): gmshfem::field::Field< std::complex<double> , gmshfem::field::Form::Form0 >(name, domain, gmshfem::field::FunctionSpaceTypeForm0(), 0, gmodel) {};
   ModelField(std::string name,const gmshfem::domain::Domain &domain, std::string gmodel, const model::Discretization& m_discret): gmshfem::field::Field< std::complex<double> , gmshfem::field::Form::Form0 >(name, domain, m_discret.functionSpaceType(), m_discret.functionSpaceDegree(), gmodel) {};
   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);}
 };
 

common/model/innerproduct/innerproduct.cpp

@@ -13,8 +13,8 @@ static const std::complex< double > im = std::complex< double >(0., 1.);
 /*
 * DifferentialInnerProductInterface
 */
-DifferentialInnerProductInterface::DifferentialInnerProductInterface(const ConfigurationInterface* const config, const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem)
-: InnerProductInterface(config, m_discret), _j("_j",_config->model_evaldomain(),_config->model_gmodel(),_m_discret), _formulation("model_formulation"),_systemIsUpToDate(false),_systemIsFactorized(false)
+DifferentialInnerProductInterface::DifferentialInnerProductInterface(bool complex, const ConfigurationInterface* const config, const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem)
+: 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"))
     {
@@ -45,15 +45,30 @@ const ModelField& DifferentialInnerProductInterface::update(const ModelFunction&
     _systemIsFactorized=true;
     _formulation.setRHSToZero();
 
+    if(!_complex)
+    {
+        for(auto it = _j.begin(); it != _j.end(); it++)
+        {
+            it->second = ((std::complex< double >)std::real(it->second));
+        }
+    }
+
     return _j;
 }
 
-std::complex<double> DifferentialInnerProductInterface::product(const ModelField& m1, const ModelField& m2, bool conjugate = false)
+std::complex<double> DifferentialInnerProductInterface::product(const ModelField& m1, const ModelField& m2)
 {
-    algebra::MatrixCCS< std::complex<double> > M;
+    algebra::MatrixCRSFast< std::complex<double> > M;
     algebra::Vector< std::complex<double> > m1v, m2v;
     m1.getVector(m1v);
     m2.getVector(m2v);
     _formulation.getLHS(M);
+    if(_complex)
+    {
+        for(unsigned int n = 0; n < m1v.size(); n++)
+        {
+            m1v[n] = std::conj(m1v[n]);
+        }
+    }
     return bilinear(m1v, M, m2v);
 }

common/model/innerproduct/innerproduct.h

@@ -20,17 +20,18 @@ class ModelUpdater;
 class InnerProductInterface
 {
 protected:
+    const bool _complex;
     const ConfigurationInterface* const _config;
     const model::Discretization _m_discret;
 public:
-    InnerProductInterface(const ConfigurationInterface* const config, const model::Discretization& m_discret): _config(config), _m_discret(m_discret) {};
+    InnerProductInterface(bool complex, const ConfigurationInterface* const config, const model::Discretization& m_discret): _complex(complex), _config(config), _m_discret(m_discret) {};
     virtual ~InnerProductInterface() = default;
 
     virtual void link(ModelUpdater* const mu, SensitivityUpdater* const su) {};
     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, bool conjugate = false) = 0;
+    virtual std::complex<double> product(const ModelField& m1, const ModelField& m2) = 0;
 
     virtual void modelIsObsolete() {};
 };
@@ -48,11 +49,11 @@ protected:
     bool _systemIsUpToDate;
     bool _systemIsFactorized;
 public:
-    DifferentialInnerProductInterface(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);
     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, bool conjugate);
+    virtual std::complex<double> product(const ModelField& m1, const ModelField& m2);
 protected:
     virtual void setLHS() = 0;
     virtual void setRHS(const ModelFunction& bulk_sensi, const ModelFunction& boundary_sensi) = 0;

common/optimization/descentsearch.h

+#ifndef H_COMMON_OPTIMIZATION_DESCENTSEARCH
+#define H_COMMON_OPTIMIZATION_DESCENTSEARCH
+
+//GmshFWI Library
+#include "../functional.h"
+
+/*
+* class DescentSearchHistoryInterface
+*/
+class DescentSearchHistoryInterface
+{
+public:
+    virtual void write(std::string filename) const = 0;
+};
+
+/*
+* class DescentSearchInterface
+*/
+/*
+ computes the descent direction p(m) where m = functional.m()
+ stores the descent direction in functional.dm()
+*/
+class DescentSearchInterface
+{
+public:
+    virtual ~DescentSearchInterface() = default;
+
+    virtual void descentsearch(FunctionalInterface* const functional) const = 0;
+
+    virtual const DescentSearchHistoryInterface* const history() const = 0;
+};
+
+#endif // H_COMMON_OPTIMIZATION_LINESEARCH

common/optimization/linesearch.h

+#ifndef H_COMMON_OPTIMIZATION_LINESEARCH
+#define H_COMMON_OPTIMIZATION_LINESEARCH
+
+//GmshFWI Library
+#include "../functional.h"
+
+/*
+* class LineSearchHistoryInterface
+*/
+class LineSearchHistoryInterface
+{
+public:
+    virtual void write(std::string filename) const = 0;
+};
+
+/*
+* class LineSearchInterface
+*/
+/*
+ returns argmin_a J(m0+a*p)
+    where J = functional
+    m0 = functional.m()
+    p = functional.dm()
+ At the end of the linesearch functional.m() is set to m+ap
+*/
+class LineSearchInterface
+{
+public:
+    virtual ~LineSearchInterface() = default;
+
+    virtual double linesearch(FunctionalInterface* const functional) const = 0;
+
+    virtual const LineSearchHistoryInterface* const history() const = 0;
+};
+
+#endif // H_COMMON_OPTIMIZATION_LINESEARCH

common/optimization/minimumsearch.h

+#ifndef H_COMMON_OPTIMIZATION_MINIMUMSEARCH
+#define H_COMMON_OPTIMIZATION_MINIMUMSEARCH
+
+//GmshFWI Library
+#include "../functional.h"
+#include "../model/element.h"
+#include "linesearch.h"
+#include "descentsearch.h"
+
+/*
+* class MinimumSearchHistoryInterface
+*/
+class MinimumSearchHistoryInterface
+{
+public:
+    virtual void write(std::string filename) const = 0;
+};
+
+/*
+* class MinimumSearchInterface
+*/
+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 const MinimumSearchHistoryInterface* const history() const = 0;
+};
+
+#endif // H_COMMON_OPTIMIZATION_MINIMUMSEARCH

gradient.cpp

@@ -55,7 +55,13 @@ int gradient(const GmshFem& gmshFem)
     ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(&d0,gmshFem);
 
     model::Discretization m_discret(gmshFem);
-    sobolev::InnerProduct innerproduct(configuration,m_discret,gmshFem);
+    unsigned int complex_buffer = 0;
+    if(!gmshFem.userDefinedParameter(complex_buffer, "complex"))
+    {
+        throw Exception("Model scalar type (complex) could not be found.");
+    }
+    bool complex = ((bool) complex_buffer);
+    sobolev::InnerProduct innerproduct(complex,configuration,m_discret,gmshFem);
 
     StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,&innerproduct,equation,objective);
     FunctionalInterface* const functional = statefunctional;

input/paper1/gradient_highly.txt

@@ -19,6 +19,9 @@ N = 2
 #Data
 data = filled_inclusion_data
 #
+#Inner product
+complex = 0
+#
 #Discretization
 h = 0.25
 model_FunctionSpaceType = HierarchicalH1

specific/model/innerproduct/newInnerProduct.h

@@ -11,14 +11,14 @@
 //GmshFWI Library
 #include "sobolev.h"
 
-InnerProductInterface* newInnerProduct(const ConfigurationInterface* const config,  const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem)
+InnerProductInterface* newInnerProduct(bool complex, const ConfigurationInterface* const config,  const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem)
 {
     std::string innerproduct;
     if(!gmshFem.userDefinedParameter(innerproduct, "innerproduct"))
     {
       throw gmshfem::common::Exception("Inner product type could not be found.");
     }
-    if(innerproduct=="sobolev"){return new sobolev::InnerProduct(config,m_discret,gmshFem);}
+    if(innerproduct=="sobolev"){return new sobolev::InnerProduct(complex,config,m_discret,gmshFem);}
     else
     {
         throw gmshfem::common::Exception("Inner product " + innerproduct + " is not valid.");

specific/model/innerproduct/sobolev.h

@@ -23,7 +23,7 @@ namespace sobolev
         double _bulkWeight;
         double _boundaryWeight;
     public:
-        InnerProduct(const ConfigurationInterface* const config, const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem) : DifferentialInnerProductInterface(config,m_discret,gmshFem), _mu(nullptr), _su(nullptr), _diag_preconditioner(false), _alpha1(0.), _bulkWeight(1.0), _boundaryWeight(0.) {};
+        InnerProduct(bool complex, const ConfigurationInterface* const config, const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem) : DifferentialInnerProductInterface(complex, config,m_discret,gmshFem), _mu(nullptr), _su(nullptr), _diag_preconditioner(false), _alpha1(0.), _bulkWeight(1.0), _boundaryWeight(0.) {};
         virtual void link(ModelUpdater* const mu, SensitivityUpdater* const mpu) override;
         virtual void unlink() override;
         virtual void modelIsObsolete() override;

specific/optimization/minimumsearch/classic.cpp

+//GmshFEM Library
+#include "Exception.h"
+//#include "Function.h" // to be removed
+
+//GmshFWI Library
+#include "classic.h"
+
+using namespace gmshfem;
+using namespace gmshfem::common;
+//using namespace gmshfem::function; // to be removed
+
+namespace classic
+{
+    /*
+    * class MinimumSearchHistoryLine2
+    */
+    void MinimumSearchHistoryLine2::write(CSVio& file, bool write_header) const
+    {
+        if(write_header)
+        {
+            file << "Performance" << "Gradient norm squared" << csv::endl;
+        }
+        else
+        {
+            file << performance << gradientNorm2 << csv::endl;
+        }
+    }
+
+    /*
+    * class MinimumSearchHistoryLine
+    */
+    void MinimumSearchHistoryLine::write(CSVio& file, bool write_header) const
+    {
+        if(write_header)
+        {
+            file << "Success" << "Iteration number" << csv::endl;
+            file << "History" << csv::endl;
+        }
+        else
+        {
+            file << success << iterations << csv::endl;
+            if(!history.empty())
+            {
+                history.begin()->write(file,true);
+                for (auto it = history.begin(); it != history.end(); it++)
+                {
+                    it->write(file);
+                }
+            }
+        }
+    }
+
+    /*
+    * class MinimumSearchHistory
+    */
+    void MinimumSearchHistory::write(std::string filename) const
+    {
+        CSVio file(filename, ';', common::OpeningMode::NewFile);
+        if(this->empty())
+        {
+            file << "Minimum search history is empty." << csv::endl;
+        }
+        else
+        {
+            this->begin()->write(file,true);
+            for (auto it = this->begin(); it != this->end(); it++)
+            {
+                file << "--" << csv::endl;
+                it->write(file);
+            }
+        }
+    }
+
+    /*
+    * class MinimumSearch
+    */
+    MinimumSearch::MinimumSearch(const GmshFem& gmshFem) : _history( new MinimumSearchHistory() )
+    {
+        if(!gmshFem.userDefinedParameter(_maxIteration, "maxIteration"))
+        {
+            throw Exception("Maximum number of iterations (maxIteration) could not be found.");
+        }
+    }
+
+    void MinimumSearch::minimize(ModelField* const m, FunctionalInterface* const functional, const DescentSearchInterface* const descentsearch, const LineSearchInterface* const linesearch) const
+    {
+        std::vector<MinimumSearchHistoryLine2> historyline;
+
+        double jn = 0.;
+        double djn = 0.;
+        unsigned int n = 0;
+        bool success = false;
+        while(true)
+        {
+            msg::print << "--- Iteration #"+ std::to_string(n)+" ----" << msg::endl;
+            jn = functional->performance(*m);
+            msg::print << "\t" << "performance = " <<  jn  << msg::endl;
+            functional->gradient().write("gradient"+std::to_string(n));
+            djn = functional->directional1(functional->gradient());
+            msg::print << "\t" << "gradient norm 2 = " <<  djn  << msg::endl;
+            historyline.emplace_back(jn,djn);
+
+            if(n >= _maxIteration){break;} // ou critère sur djn ou critère sur jn
+            else
+            {
+                //descent_direction();
+                //linesearch();
+                //m = m+a*p
+            }
+            n++;
+        }
+
+        _history->emplace_back(success,n,historyline);
+    }
+}

specific/optimization/minimumsearch/classic.h

+#ifndef H_SPECIFIC_OPTIMIZATION_MINIMUMSEARCH_CLASSIC
+#define H_SPECIFIC_OPTIMIZATION_MINIMUMSEARCH_CLASSIC
+
+//GmshFEM Library
+#include "GmshFem.h"
+#include "CSVio.h"
+
+//GmshFWI Library
+#include "../../../common/optimization/minimumsearch.h"
+
+namespace classic
+{
+    /*
+    * class MinimumSearchHistoryLine2
+    */
+    class MinimumSearchHistoryLine2
+    {
+    public:
+        const double performance;
+        const double gradientNorm2;
+
+        MinimumSearchHistoryLine2(double performance_init, double gradientNorm2_init) : performance(performance_init), gradientNorm2(gradientNorm2_init) {};
+
+        void write(gmshfem::common::CSVio& file, bool write_header=false) const;
+    };
+
+    /*
+    * class MinimumSearchHistoryLine
+    */
+    class MinimumSearchHistoryLine
+    {
+    public:
+        const bool success;
+        const unsigned int iterations;
+        const std::vector<MinimumSearchHistoryLine2> history;
+
+        MinimumSearchHistoryLine(bool success_init, unsigned int iterations_init, const std::vector<MinimumSearchHistoryLine2>& history_init) : success(success_init), iterations(iterations_init), history(history_init) {};
+
+        void write(gmshfem::common::CSVio& file, bool write_header=false) const;
+    };
+
+    /*
+    * class MinimumSearchHistory
+    */
+    class MinimumSearchHistory final: public MinimumSearchHistoryInterface, public std::vector<MinimumSearchHistoryLine>
+    {
+    public:
+        virtual void write(std::string filename) const;
+    };
+
+    /*
+    * class MinimumSearch
+    */
+    class MinimumSearch final: public MinimumSearchInterface
+    {
+    private:
+        unsigned int _maxIteration;
+        MinimumSearchHistory* const _history;
+    public:
+        MinimumSearch(const gmshfem::common::GmshFem& gmshFem);
+        ~MinimumSearch() {delete _history;};
+
+        virtual void minimize(ModelField* const m, FunctionalInterface* const functional, const DescentSearchInterface* const descentsearch, const LineSearchInterface* const linesearch) const;
+
+        virtual const MinimumSearchHistoryInterface* const history() const {return _history;};
+    };
+};
+
+#endif // H_SPECIFIC_OPTIMIZATION_MINIMUMSEARCH_CLASSIC

specific/optimization/minimumsearch/newMinimumSearch.h

+#ifndef H_SPECIFIC_WAVE_NEWMINIMUMSEARCH
+#define H_SPECIFIC_WAVE_NEWMINIMUMSEARCH
+
+//Standard Library
+#include <string>
+
+//GmshFem Library
+#include "GmshFem.h"
+#include "Exception.h"
+
+//FWI Library
+#include "classic.h"
+
+MinimumSearchInterface* newMinimumSearch(const gmshfem::common::GmshFem& gmshFem)
+{
+    std::string minimumsearch;
+    if(!gmshFem.userDefinedParameter(minimumsearch, "minimumsearch"))
+    {
+      throw gmshfem::common::Exception("Minimum search type could not be found.");
+    }
+    if(minimumsearch=="classic"){return new classic::MinimumSearch(gmshFem);}
+    else
+    {
+        throw gmshfem::common::Exception("Minimum search " + minimumsearch + " is not valid.");
+        return nullptr;
+    }
+};
+
+#endif //H_SPECIFIC_WAVE_NEWMINIMUMSEARCH