Solve memory leak (m0 velocity soil)

common/configuration.h

@@ -60,6 +60,7 @@ public:
     virtual std::string data_gmodel() const = 0;
 
     virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const = 0;
+    virtual std::array<double,2> index_to_data_coordinate(unsigned int s, unsigned int r) const = 0;
 
     gmshfem::domain::Domain wave_omega(Support support) const {return _wave_omega[support];};
     gmshfem::domain::Domain model_unknown(Support support) const {return _model_unknown[support];};
@@ -90,7 +91,7 @@ public:
     std::vector<unsigned int> receiverPerShot() const;
     ModelFunction m0() const
     {
-      return _m0.copy();
+      return _m0;
     };
 
     virtual double m_reference() const = 0;

common/data/element.cpp

@@ -7,7 +7,7 @@ using namespace gmshfem::common;
 
 //FWI Library
 #include "element.h"
-#include "../../specific/wavetodata.h"
+#include "../../specific/todata.h"
 
 /*
 * Class Data
@@ -121,6 +121,11 @@ void Data<T_Physic>::value(const WaveMultiField<T_Physic>& w)
     *this = wave_to_data(_config,w);
 }
 
+template<Physic T_Physic>
+void Data<T_Physic>::value(const DataField<T_Physic>& d)
+{
+    *this = data_to_data(_config,d);
+}
 
 template<Physic T_Physic>
 std::vector<unsigned int> Data<T_Physic>::receiverPerShot() const

common/data/element.h

@@ -27,6 +27,7 @@ public:
     PointData<T_Physic> value(double  xs, double xr) const;
     void value(unsigned int s,unsigned int r,const PointData<T_Physic>& value);
     void value(const WaveMultiField<T_Physic>& w);
+    void value(const DataField<T_Physic>& d);
     unsigned int ns() const {return _value.size();};
     unsigned int nr(unsigned int shot) const {return _value[shot].size();};
     std::vector<unsigned int> receiverPerShot() const;

common/model/macro.cpp

+//GmshFEM Library
+#include "Formulation.h"
+#include "Function.h"
+
+//GmshFWI Library
+#include "macro.h"
+
+using namespace gmshfem::problem;
+using namespace gmshfem::equation;
+using namespace gmshfem::function;
+
+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;
+}
+
+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;
+}
+
+ModelField laplace_filter(const ModelFunction& m0, double scale, const ConfigurationInterface* const config, const model::Discretization& m_discret, std::string integrationType)
+{
+    double alpha1 = (scale / (2 * M_PI)) * (scale / (2 * M_PI));
+    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 Laplace filter");
+    formulation.galerkin(dof(m), tf(m),config->model_unknown(Support::BLK), integrationType);
+    formulation.galerkin(alpha1*grad(dof(m)), grad(tf(m)),config->model_unknown(Support::BLK), integrationType);
+    formulation.galerkin(-m0,tf(m), config->model_unknown(Support::BLK), integrationType);
+    //formulation.galerkin(-1.,tf(m), config->model_unknown(Support::BLK), integrationType);
+    formulation.pre();
+    formulation.assemble();
+    formulation.solve();
+    return m;
+}

common/model/macro.h

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

common/statefunctional.cpp

@@ -159,9 +159,9 @@ const ModelField& StateFunctional<T_Physic>::gradient()
     return _mu.get(ToBeUpdated).state(Type::A).field();
 }
 
-/* directional */
+/* sensitivity */
 template<Physic T_Physic>
-double StateFunctional<T_Physic>::directional(Order order, Support support, const ModelFunction &dm2)
+const ModelFunction& StateFunctional<T_Physic>::sensitivity(Order order, Support support)
 {
     std::array<bool,4> ModelToBeUpToDate = {false,false,false,false};
     switch (order)
@@ -171,10 +171,17 @@ double StateFunctional<T_Physic>::directional(Order order, Support support, cons
         case Order::FST: case Order::DIAG:
             ModelToBeUpToDate[Type::F] = true;
     }
+    return _su.get(order,support,_mu.get(ModelToBeUpToDate));
+}
+
+/* directional */
+template<Physic T_Physic>
+double StateFunctional<T_Physic>::directional(Order order, Support support, const ModelFunction &dm2)
+{
     return std::real(
         _equation->directional(
             support,
-            _su.get(order,support,_mu.get(ModelToBeUpToDate)),
+            sensitivity(order,support),
             dm2
         )
     );

common/statefunctional.h

@@ -55,6 +55,8 @@ public:
     /* gradient */
     virtual const ModelField& gradient();
 
+    /* sensitivity */
+    const ModelFunction& sensitivity(Order order, Support support);
     /* directional */
     double directional(Order order , Support support, const ModelFunction &dm);
 

common/wave/preconditioner.oldcpp

-//Standard library
-
-//GmshFem library
-#include "Exception.h"
-
-//GmshFWI Library
-#include "preconditioner.h"
-
-using namespace gmshfem;
-using namespace gmshfem::common;
-
-/*
-* PreconditionerState
-*/
-const ModelFunction& ModelPreconditionerState::state() const
-{
-    if(!_isUpToDate)
-    {
-        throw Exception("Model preconditioner is not up to date while required.");
-    }
-    return _state;
-}
-
-/*
-* ModelPreconditionerUpdater
-*/
-void ModelPreconditionerUpdater::update(const ModelStateEvaluator& ms)
-{
-    if(_ps._isUpToDate){return;}
-    _ps._state = update_preconditioner(ms);
-    _ps._isUpToDate = true;
-}
-
-
-const ModelFunction& ModelPreconditionerUpdater::get_preconditioner(const ModelStateEvaluator& ms)
-{
-    update(ms);
-    return _ps.state();
-}
-
-void ModelPreconditionerUpdater::isObsolete(bool NoMoreUpToDate)
-{
-    _ps._isUpToDate = !NoMoreUpToDate;
-}

common/wave/preconditioner.oldh

-#ifndef H_COMMON_WAVE_PRECONDITIONER
-#define H_COMMON_WAVE_PRECONDITIONER
-
-//FWI Library
-#include "../model/state.h"
-
-/*
-* ModelPreconditionerState
-*/
-class ModelPreconditionerState
-{
-private:
-    ModelFunction _state;
-    bool _isUpToDate;
-public:
-    const ModelFunction& state() const;
-
-friend class ModelPreconditionerUpdater;
-};
-
-/*
-* ModelPreconditionerUpdater
-*/
-class ModelPreconditionerUpdater
-{
-private:
-    ModelPreconditionerState _ps;
-    void update(const ModelStateEvaluator& ms);
-    virtual ModelFunction update_preconditioner(const ModelStateEvaluator& ms) = 0;
-public:
-    const ModelFunction& get_preconditioner(const ModelStateEvaluator& ms);
-    void isObsolete(bool NoMoreUpToDate);
-    virtual bool isModelDependent() {return true;};
-};
-
-#endif // H_COMMON_WAVE_PRECONDITIONER

inversion.cpp

@@ -10,8 +10,7 @@ using namespace gmshfem::common;
 #include "specific/wave/equation/newEquation.h"
 #include "specific/data/objective/newObjective.h"
 #include "specific/model/innerproduct/newInnerProduct.h"
-#include "common/initialization.h"
-#include "common/rediscretize.h"
+#include "common/model/macro.h"
 #include "common/statefunctional.h"
 #include "specific/optimization/minimumsearch/newMinimumSearch.h"
 #include "specific/optimization/descentsearch/newDescentSearch.h"

specific/configuration/inclusion.cpp

@@ -260,7 +260,7 @@ namespace inclusion
         ScalarPiecewiseFunction< std::complex< double > > m0;
         m0.addFunction(_mb,_background[Support::BLK] | _background[Support::BND] | _points);
         m0.addFunction(_mc,_inclusion[Support::BLK] | _inclusion[Support::BND]);
-        _m0 = m0.copy();
+        _m0 = m0;
     }
 
     void Configuration::mesh() const
@@ -514,6 +514,11 @@ namespace inclusion
         throw Exception("Data coordinate to index not yet implemented for inclusion.");
     }
 
+    std::array<double,2> Configuration::index_to_data_coordinate(unsigned int s, unsigned int r) const
+    {
+        throw Exception("Index to data coordinates not yet implemented for inclusion.");
+    }
+
     template<Physic T_Physic>
     ModelFunction green0_preconditioner(const WaveMultiField<T_Physic>& g, const Configuration* const config)
     {

specific/configuration/inclusion.h

@@ -87,6 +87,7 @@ namespace inclusion
         virtual void mesh() const;
 
         virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const;
+        virtual std::array<double,2> index_to_data_coordinate(unsigned int s, unsigned int r) const;
 
         virtual double area() const {return _L * _H; };
         virtual double array() const {return std::max(_He,_Le);};

specific/configuration/soil.cpp

@@ -244,7 +244,7 @@ namespace soil
           m0.addFunction(1. / pow(vy,2),_soil[Support::BLK] | _soil[Support::BND]);
           m0.addFunction(_m_air, _air[Support::BLK] | _air[Support::BND]);
         }
-        _m0 = m0.copy();
+        _m0 = m0;
     }
 
     void Configuration::mesh() const
@@ -425,6 +425,16 @@ namespace soil
         return index;
     }
 
+    std::array<double,2> Configuration::index_to_data_coordinate(unsigned int s, unsigned int r) const
+    {
+        std::array<double,2> coord;
+        double x0 = (_L-array())/2.;
+        double step = array()/(ner()-1);
+        coord[0] = x0 + s*step;
+        coord[1] = x0 + r*step;
+        return coord;
+    }
+
     template<Physic T_Physic>
     ModelFunction green0_preconditioner(const WaveMultiField<T_Physic>& g, const Configuration* const config)
     {

specific/configuration/soil.h

@@ -46,6 +46,7 @@ namespace soil
     unsigned int ner() const {return _ner;};
 
     virtual std::array<unsigned int,2> data_coordinate_to_index(double xs, double xr) const;
+    virtual std::array<double,2> index_to_data_coordinate(unsigned int s, unsigned int r) const;
 
     double H() const {return _H;};
     virtual double area() const {return _L * _H;};

specific/data/element.cpp

@@ -42,7 +42,7 @@ void DataField<Physic::acoustic>::forward_interpolate(const Data<Physic::acousti
         case 1:
             for(auto it = this->begin(); it != this->end(); it++)
             {
-                const dofs::Dof< scalar::Precision< std::complex<double> > > *dof = it->first;
+                const dofs::Dof *dof = it->first;
                 it->second = data.value(dof->x(),dof->y());
             }
             break;

specific/wavetodata.cpp → specific/todata.cpp

@@ -2,8 +2,8 @@
 #include "Integrate.h"
 #include "Function.h"
 
-//FWI Library
-#include "wavetodata.h"
+//GmshFWI Library
+#include "todata.h"
 
 using namespace gmshfem::post;
 using namespace gmshfem::function;
@@ -22,3 +22,21 @@ Data<Physic::acoustic> wave_to_data(const ConfigurationInterface* const config,
         }
         return v;
 }
+
+template<>
+Data<Physic::acoustic> data_to_data(const ConfigurationInterface* const config, const DataField<Physic::acoustic>& d)
+{
+        Data<Physic::acoustic> v(config);
+        for (unsigned int s = 0; s < config->ns(); s++)
+        {
+            for (unsigned int r = 0; r < config->nr(s); r++)
+            {
+                    std::array<double,2> coord = config->index_to_data_coordinate(s,r);
+                    double xs = coord[0];
+                    double xr = coord[1];
+                    PointData<Physic::acoustic> buffer(evaluate(d,xs,xr,0.));
+                    v.value(s,r,buffer);
+            }
+        }
+        return v;
+}

specific/wavetodata.h → specific/todata.h

-#ifndef H_SPECIFIC_RESTRICTION
-#define H_SPECIFIC_RESTRICTION
+#ifndef H_SPECIFIC_TODATA
+#define H_SPECIFIC_TODATA
 
 #include "../common/data/element.h"
 #include "../common/wave/element.h"
@@ -7,4 +7,7 @@
 template<Physic T_Physic>
 Data<T_Physic> wave_to_data(const ConfigurationInterface* const config, const WaveMultiField<T_Physic>& w);
 
-#endif // H_SPECIFIC_RESTRICTION
+template<Physic T_Physic>
+Data<T_Physic> data_to_data(const ConfigurationInterface* const config, const DataField<T_Physic>& d);
+
+#endif // H_SPECIFIC_TODATA

synthetics.cpp

@@ -46,7 +46,7 @@ int synthetics(const GmshFem& gmshFem)
 
         std::array<bool,4> needUpToDate = {true,false,false,false};
         ((du.get(needUpToDate,ms)).state(Type::F)).write(name+"_data"+suffix);
-        ((du.get(needUpToDate,ms)).state(Type::F)).write(name+"_data"+suffix,"pos");
+        ((du.get(needUpToDate,ms)).state(Type::F)).write(name+"_data_field"+suffix,"pos");
 
         unsigned int write_fields = 1;
         gmshFem.userDefinedParameter(write_fields, "write_fields");