uptodate

CMakeLists.txt

@@ -109,15 +109,15 @@ add_executable( directional directional.cpp ${EXTRA_INCS})
 target_link_libraries( directional fwi  ${EXTRA_LIBS})
 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( preconditioner preconditioner.cpp ${EXTRA_INCS})
-target_link_libraries( preconditioner fwi  ${EXTRA_LIBS})
+#add_executable( convergence convergence.cpp ${EXTRA_INCS})
+#target_link_libraries( convergence fwi  ${EXTRA_LIBS})
+#add_executable( preconditioner preconditioner.cpp ${EXTRA_INCS})
+#target_link_libraries( preconditioner fwi  ${EXTRA_LIBS})
 add_executable( ip_comparison ip_comparison.cpp ${EXTRA_INCS})
 target_link_libraries( ip_comparison fwi  ${EXTRA_LIBS})
-add_executable( ob_comparison ob_comparison.cpp ${EXTRA_INCS})
-target_link_libraries( ob_comparison fwi  ${EXTRA_LIBS})
-add_executable( multiscale multiscale.cpp ${EXTRA_INCS})
-target_link_libraries( multiscale fwi  ${EXTRA_LIBS})
+#add_executable( ob_comparison ob_comparison.cpp ${EXTRA_INCS})
+#target_link_libraries( ob_comparison fwi  ${EXTRA_LIBS})
+#add_executable( multiscale multiscale.cpp ${EXTRA_INCS})
+#target_link_libraries( multiscale fwi  ${EXTRA_LIBS})
 add_executable( inversion inversion.cpp ${EXTRA_INCS})
 target_link_libraries( inversion fwi  ${EXTRA_LIBS})

common/model/innerproduct/innerproduct.cpp

@@ -26,20 +26,20 @@ DifferentialInnerProductInterface::DifferentialInnerProductInterface(const Confi
       throw Exception("Model integration type could not be found.");
     }
     if(!(
-        gmshFem.userDefinedParameter(_integrationDegreeBlk, "integration_degree_blk")
+        gmshFem.userDefinedParameter(_integrationDegreeBlk, "innerproduct_integration_degree_blk")
         ||
-        gmshFem.userDefinedParameter(_integrationDegreeBlk, "integration_degree_blk"+suffix)
+        gmshFem.userDefinedParameter(_integrationDegreeBlk, "innerproduct_integration_degree_blk"+suffix)
     ))
     {
-      throw Exception("Model integration degree (bulk) could not be found.");
+      throw Exception("Inner product integration degree (bulk) could not be found.");
     }
     if(!(
-        gmshFem.userDefinedParameter(_integrationDegreeBnd, "integration_degree_bnd")
+        gmshFem.userDefinedParameter(_integrationDegreeBnd, "innerproduct_integration_degree_bnd")
         ||
-        gmshFem.userDefinedParameter(_integrationDegreeBnd, "integration_degree_bnd"+suffix)
+        gmshFem.userDefinedParameter(_integrationDegreeBnd, "innerproduct_integration_degree_bnd"+suffix)
     ))
     {
-      throw Exception("Model integration degree (boundary) could not be found.");
+      throw Exception("Inner product integration degree (boundary) could not be found.");
     }
 }
 

common/optimization/globalminimumsearch.h

@@ -28,7 +28,7 @@ public:
     virtual ~GlobalMinimumSearchInterface() = default;
     void name(const std::string& name){_name=name;};
 
-    virtual void operator()(ModelField* const m, const Data<T_Physic>* const d0, EquationInterface<T_Physic>* const equation) const = 0;
+    virtual void operator()(ModelField* const m, const std::vector<const Data<T_Physic>*>& d0, const std::vector<EquationInterface<T_Physic>*>& equation) const = 0;
 
     virtual const GlobalMinimumSearchHistoryInterface* const history() const = 0;
 };

common/optimization/localminimumsearch.h

@@ -64,7 +64,7 @@ public:
         if(!gmshFem.userDefinedParameter(_overallRelDecreaseThreshold, "localminimum_overallRelDecreaseThreshold"))
         {
             gmshfem::msg::warning << "Overall realtive decrease threshold (localminimum_overallDecreaseThreshold) could not be found." << gmshfem::msg::endl;
-            _overallRelDecreaseThreshold = 1.;
+            _overallRelDecreaseThreshold = 0.;
         }
         if(!gmshFem.userDefinedParameter(_meanRelDecreaseThreshold, "localminimum_meanRelDecreaseThreshold"))
         {

common/sensitivity/updater.cpp

@@ -61,4 +61,17 @@ Sensitivity SensitivityUpdaterFromEquation<T_Physic>::sensitivity(Order order, S
     return _equation->update_sensitivity(order,support,_du->get(dataNeedToBeUpToDate,ms),ms,_wu->get(waveNeedToBeUpToDate,ms));
 }
 
+/*
+* SensitivityUpdaterFromVector
+*/
+Sensitivity SensitivityUpdaterFromVector::sensitivity(Order order, Support support, const ModelStateEvaluator& ms)
+{
+    Sensitivity sum = _sue[0]->get(order,support,ms);
+    for (unsigned int f = 1; f < _nf; f++)
+    {
+        sum = sum + _sue[f]->get(order,support,ms);
+    }
+    return sum;
+}
+
 template class SensitivityUpdaterFromEquation<Physic::acoustic>;

common/sensitivity/updater.h

@@ -18,6 +18,8 @@ private:
     void update(Order order, Support support, const ModelStateEvaluator& ms);
     virtual Sensitivity sensitivity(Order order, Support support, const ModelStateEvaluator& ms) = 0;
 public:
+    virtual ~SensitivityUpdater() = default;
+
     const Sensitivity& get(Order order, Support support,const ModelStateEvaluator& ms);
     void isObsolete(std::array<bool,3> NoMoreUpToDate);
 };
@@ -38,4 +40,18 @@ private:
     virtual Sensitivity sensitivity(Order order, Support support, const ModelStateEvaluator& ms);
 };
 
+/*
+* SensitivityUpdaterFromVector
+*/
+class SensitivityUpdaterFromVector : public SensitivityUpdater
+{
+public:
+    SensitivityUpdaterFromVector(unsigned int nf, SensitivityUpdater* *const sue) : _nf(nf), _sue(sue) {};
+private:
+    unsigned int _nf;
+    SensitivityUpdater* *const _sue;
+
+    virtual Sensitivity sensitivity(Order order, Support support, const ModelStateEvaluator& ms);
+};
+
 #endif // H_COMMON_SENSITIVITY_UPDATER

common/statefunctional.cpp

@@ -2,17 +2,41 @@
 #include "statefunctional.h"
 
 template<Physic T_Physic>
-StateFunctional<T_Physic>::StateFunctional(const ConfigurationInterface* const config, InnerProductInterface* const innerproduct, EquationInterface<T_Physic>* const equation,  ObjectiveInterface<T_Physic>* const objective): _innerproduct(innerproduct), _equation(equation), _objective(objective), _mu(config,&_su,_innerproduct), _su(&_wu,&_du,_equation), _wu(config,&_du,_equation), _du(config,&_wu,_objective)
+StateFunctional<T_Physic>::StateFunctional(const ConfigurationInterface* const config, InnerProductInterface* const innerproduct, const std::vector<EquationInterface<T_Physic>*>& equation, const std::vector<ObjectiveInterface<T_Physic>*>& objective): _innerproduct(innerproduct), _equation(equation), _objective(objective), _nf(_equation.size()), _mu(config,&_su,_innerproduct), _sue(new SensitivityUpdater*[_nf]), _su(_nf,_sue), _wu( (WaveUpdater<T_Physic> *) std::malloc(_nf * sizeof(WaveUpdater<T_Physic>))), _du( (DataUpdater<T_Physic> *) std::malloc(_nf * sizeof(DataUpdater<T_Physic>)))
 {
+    for (unsigned int f = 0; f < _nf; f++)
+    {
+        new(&_wu[f]) WaveUpdater<T_Physic>(config,&_du[f],_equation[f]);
+        if(_objective.size()==_nf)
+        {
+            new(&_du[f]) DataUpdater<T_Physic>(config,&_wu[f],_objective[f]);
+            _objective[f]->link(&_mu,&_du[f]);
+        }
+        else if(_objective.size()==0)
+        {
+            new(&_du[f]) DataUpdater<T_Physic>(config,&_wu[f],nullptr);
+        }
+        else
+        {
+            throw gmshfem::common::Exception("Objectives vector has not the same size than equations vector");
+        }
+        _sue[f] = new SensitivityUpdaterFromEquation<T_Physic>(&_wu[f],&_du[f],_equation[f]);
+    }
     if(_innerproduct!=nullptr){_innerproduct->link(&_mu,&_su);}
-    if(_objective!=nullptr){_objective->link(&_mu,&_du);}
 }
 
 template<Physic T_Physic>
 StateFunctional<T_Physic>::~StateFunctional()
 {
     if(_innerproduct!=nullptr){_innerproduct->unlink();}
-    if(_objective!=nullptr){_objective->unlink();}
+    for (unsigned int f = 0; f < _nf; f++)
+    {
+        if(_objective.size()==_nf){_objective[f]->unlink();}
+        delete _sue[f];
+    }
+    delete[] _sue;
+    std::free(_wu);
+    std::free(_du);
 }
 
 template<Physic T_Physic>
@@ -37,7 +61,11 @@ void StateFunctional<T_Physic>::write_data(Type type, std::string name, std::str
         case Type::F: case Type::A: case Type::D:
             ModelNeedToBeUpToDate[Type::F] = true;
     }
-    _du.get(DataNeedToBeUpToDate,_mu.get(ModelNeedToBeUpToDate)).write(type,name,fmt);
+    for (unsigned int f = 0; f < _nf; f++)
+    {
+        std::string suffix_f = "f" + std::to_string(f);
+        _du[f].get(DataNeedToBeUpToDate,_mu.get(ModelNeedToBeUpToDate)).write(type,name+suffix_f,fmt);
+    }
 }
 
 template<Physic T_Physic>
@@ -58,7 +86,11 @@ void StateFunctional<T_Physic>::write_wave(Type type, std::string name)
             throw gmshfem::common::Exception("Impossible to write wave other than F,PF,A,PA");
 
     }
-    _wu.get(WaveNeedToBeUpToDate,_mu.get(ModelNeedToBeUpToDate)).write(type,name);
+    for (unsigned int f = 0; f < _nf; f++)
+    {
+        std::string suffix_f = "f" + std::to_string(f);
+        _wu[f].get(WaveNeedToBeUpToDate,_mu.get(ModelNeedToBeUpToDate)).write(type,name+suffix_f);
+    }
 }
 
 template<Physic T_Physic>
@@ -72,57 +104,69 @@ template<Physic T_Physic>
 void StateFunctional<T_Physic>::objectiveIsObsolete()
 {
     std::array<bool,5> dataNoMoreUpToDate = {false, false, true, true, true};
-    _du.isObsolete(dataNoMoreUpToDate);
-
     std::array<bool,4> waveNoMoreUpToDate = {false, false, true, true};
-    _wu.isObsolete(waveNoMoreUpToDate);
+    std::array<bool,3> sensitivityNoMoreUpToDate = {true, true, true};
+    for (unsigned int f = 0; f < _nf; f++)
+    {
+        _du[f].isObsolete(dataNoMoreUpToDate);
+        _wu[f].isObsolete(waveNoMoreUpToDate);
+        _sue[f]->isObsolete(sensitivityNoMoreUpToDate);
+    }
+    _su.isObsolete(sensitivityNoMoreUpToDate);
 
     std::array<bool,4> modelNoMoreUpToDate = {false, false, true, true};
     _mu.isObsolete(modelNoMoreUpToDate);
-
-    std::array<bool,3> sensitivityNoMoreUpToDate = {true, true, true};
-    _su.isObsolete(sensitivityNoMoreUpToDate);
 }
 
 template<Physic T_Physic>
 void StateFunctional<T_Physic>::modelIsObsolete()
 {
     std::array<bool,5> dataNoMoreUpToDate = {true, true, true, true, true};
+    std::array<bool,4> waveNoMoreUpToDate = {true, true, true, true};
+    std::array<bool,3> sensitivityNoMoreUpToDate = {true, true, true};
+    bool sensitivityNoMoreUpToDateOverallDiag = false;
+
+    for (unsigned int f = 0; f < _nf; f++)
+    {
         dataNoMoreUpToDate[Type::D] = false;
-    if(_objective != nullptr){dataNoMoreUpToDate[Type::D] = _objective->modelIsObsolete();}
-    _du.isObsolete(dataNoMoreUpToDate);
+        if(_objective.size()==_nf){dataNoMoreUpToDate[Type::D] = _objective[f]->modelIsObsolete();}
+        _du[f].isObsolete(dataNoMoreUpToDate);
 
+        _wu[f].isObsolete(waveNoMoreUpToDate);
 
-    std::array<bool,4> waveNoMoreUpToDate = {true, true, true, true};
-    _wu.isObsolete(waveNoMoreUpToDate);
+        sensitivityNoMoreUpToDate[Order::DIAG]=false;
+        sensitivityNoMoreUpToDate[Order::DIAG] = _equation[f]->modelIsObsolete();
+        if(sensitivityNoMoreUpToDate[Order::DIAG]){sensitivityNoMoreUpToDateOverallDiag=true;}
+        _sue[f]->isObsolete(sensitivityNoMoreUpToDate);
+    }
+    sensitivityNoMoreUpToDate[Order::DIAG]=sensitivityNoMoreUpToDateOverallDiag;
+    _su.isObsolete(sensitivityNoMoreUpToDate);
 
     std::array<bool,4> modelNoMoreUpToDate = {true, true, true, true};
     modelNoMoreUpToDate[Type::PF]=false;
     if(_innerproduct != nullptr){_innerproduct->modelIsObsolete();}
     _mu.isObsolete(modelNoMoreUpToDate);
-
-    std::array<bool,3> sensitivityNoMoreUpToDate = {true, true, true};
-    sensitivityNoMoreUpToDate[Order::DIAG]=false;
-    if(_equation != nullptr){sensitivityNoMoreUpToDate[Order::DIAG] = _equation->modelIsObsolete();}
-    _su.isObsolete(sensitivityNoMoreUpToDate);
 }
 
 template<Physic T_Physic>
 void StateFunctional<T_Physic>::modelPerturbationIsObsolete()
 {
     std::array<bool,5> dataNoMoreUpToDate = {false, true, false, true, false};
-    _du.isObsolete(dataNoMoreUpToDate);
-
-
     std::array<bool,4> waveNoMoreUpToDate = {false, true, false, true};
-    _wu.isObsolete(waveNoMoreUpToDate);
-    if(_equation != nullptr){_equation->modelPerturbationIsObsolete();}
+    std::array<bool,3> sensitivityNoMoreUpToDate = {false, true, false};
+
+    for (unsigned int f = 0; f < _nf; f++)
+    {
+        _du[f].isObsolete(dataNoMoreUpToDate);
+        _wu[f].isObsolete(waveNoMoreUpToDate);
+        _equation[f]->modelPerturbationIsObsolete();
+        _sue[f]->isObsolete(sensitivityNoMoreUpToDate);
+    }
+    _su.isObsolete(sensitivityNoMoreUpToDate);
 
     std::array<bool,4> modelNoMoreUpToDate = {false, true, false, true};
     _mu.isObsolete(modelNoMoreUpToDate);
 
-    std::array<bool,3> sensitivityNoMoreUpToDate = {false, true, false};
-    _su.isObsolete(sensitivityNoMoreUpToDate);
 }
 
 template<Physic T_Physic>
@@ -175,9 +219,18 @@ double StateFunctional<T_Physic>::performance()
 {
     std::array<bool,5> DataToBeUpdated = {true,false,false,false,false};
     std::array<bool,4> ModelToBeUpdated = {true,false,false,false};
-    return _objective->performance(
-        _du.get(DataToBeUpdated,_mu.get(ModelToBeUpdated)).state(Type::F)
+    double sum = 0.;
+    //gmshfem::msg::print << "performance [f] = ";
+    for (unsigned int f = 0; f < _nf; f++)
+    {
+        double pf = _objective[f]->performance(
+            _du[f].get(DataToBeUpdated,_mu.get(ModelToBeUpdated)).state(Type::F)
         );
+        //gmshfem::msg::print << pf << ";";
+        sum += pf;
+    }
+    //gmshfem::msg::print << gmshfem::msg::endl;
+    return sum;
 }
 
 /*
@@ -206,18 +259,37 @@ const Sensitivity& StateFunctional<T_Physic>::sensitivity(Order order, Support s
     return _su.get(order,support,_mu.get(ModelToBeUpToDate));
 }
 
+template<Physic T_Physic>
+const Sensitivity& StateFunctional<T_Physic>::sensitivity(unsigned int f, Order order, Support support)
+{
+    std::array<bool,4> ModelToBeUpToDate = {false,false,false,false};
+    switch (order)
+    {
+        case Order::SCD:
+            ModelToBeUpToDate[Type::PF] = true;
+        case Order::FST: case Order::DIAG:
+            ModelToBeUpToDate[Type::F] = true;
+    }
+    return _sue[f]->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(
-            sensitivity(order,support),
+    double sum = 0.;
+    for (unsigned int f = 0; f < _nf; f++)
+    {
+        sum += std::real(
+            _equation[f]->directional(
+                sensitivity(f,order,support),
                 dm2,
                 support
             )
         );
     }
+    return sum;
+}
 
 /* directional1 */
 template<Physic T_Physic>

common/statefunctional.h

@@ -15,18 +15,22 @@ class StateFunctional: public FunctionalInterface
 {
 private:
     InnerProductInterface* const _innerproduct;
-    EquationInterface<T_Physic>* const _equation;
-    ObjectiveInterface<T_Physic>* const _objective;
+    const std::vector<EquationInterface<T_Physic>*> _equation;
+    const std::vector<ObjectiveInterface<T_Physic>*> _objective;
+    const unsigned int _nf;
 
     ModelUpdater _mu;
-    SensitivityUpdaterFromEquation<T_Physic> _su;
-    WaveUpdater<T_Physic> _wu;
-    DataUpdater<T_Physic> _du;
+
+    SensitivityUpdater* *const _sue;
+    SensitivityUpdaterFromVector _su;
+
+    WaveUpdater<T_Physic> *const _wu;
+    DataUpdater<T_Physic> *const _du;
 
     void modelIsObsolete();
     void modelPerturbationIsObsolete();
 public:
-    StateFunctional(const ConfigurationInterface* const config, InnerProductInterface* const innerproduct, EquationInterface<T_Physic>* const equation, ObjectiveInterface<T_Physic>* const _objective);
+    StateFunctional(const ConfigurationInterface* const config, InnerProductInterface* const innerproduct, const std::vector<EquationInterface<T_Physic>*>& equation, const std::vector<ObjectiveInterface<T_Physic>*>& objective);
 
     ~StateFunctional();
 
@@ -58,6 +62,7 @@ public:
 
     /* sensitivity */
     const Sensitivity& sensitivity(Order order, Support support);
+    const Sensitivity& sensitivity(unsigned int f, Order order, Support support);
 
     /* directional */
     double directional(Order order , Support support, const ModelFunction &dm);

common/wave/element.cpp

@@ -18,10 +18,11 @@ void WaveMultiField<T_Physic>::initializeField(const WaveField<T_Physic>& field)
 template<Physic T_Physic>
 const WaveField<T_Physic>& WaveMultiField<T_Physic>::operator[](const unsigned int index) const
 {
-    if(index!=_index)
+    if(index!=_index || !_isUpToDate)
     {
         _field.setAllVector(_dofValues[index], gmshfem::dofs::RawOrder::Hash);
         _index=index;
+        _isUpToDate=true;
     }
     return _field;
 }
@@ -60,6 +61,7 @@ void WaveAuxilaryField<T_Physic>::assignValues(const std::vector< std::complex<d
     for(unsigned int i = 0; i < values.size(); ++i) {
       _wmf->_dofValues[_index][i] = values[_systemValuesToDofValues[i]];
     }
+    if(_index == _wmf->_index){_wmf->_isUpToDate=false;}
 }
 
 template class WaveMultiField<Physic::acoustic>;

common/wave/element.h

@@ -18,10 +18,11 @@ class WaveMultiField
 private:
     mutable WaveField<T_Physic> _field;
     bool _isInitialized;
+    mutable bool _isUpToDate;
     std::vector< gmshfem::algebra::Vector< std::complex< double > > > _dofValues;
     mutable unsigned int _index;
 public:
-    WaveMultiField(unsigned int n, const std::string& name = "") : _field(), _isInitialized(false), _dofValues(n), _index(n+1) {_field.name(name);};
+    WaveMultiField(unsigned int n, const std::string& name = "") : _field(), _isInitialized(false), _isUpToDate(false), _dofValues(n), _index(n+1) {_field.name(name);};
 
     void initializeField(const WaveField<T_Physic>& field);
     bool isInitialized() const {return _isInitialized;};

common/wave/equation/equation.cpp

 //GmshFem Library
-#include "FieldRoutines.h"
+//#include "FieldRoutines.h"
 
 //FWI Library
 #include "equation.h"
+#include "../../../specific/wave/equation/newEquation.h"
 
 using namespace gmshfem;
 using namespace gmshfem::common;
@@ -21,20 +22,20 @@ EquationInterface<T_Physic>::EquationInterface(const ConfigurationInterface* con
       throw Exception("Integration type could not be found.");
     }
     if(!(
-        gmshFem.userDefinedParameter(_integrationDegreeBlk, "integration_degree_blk")
+        gmshFem.userDefinedParameter(_integrationDegreeBlk, "equation_integration_degree_blk")
         ||
-        gmshFem.userDefinedParameter(_integrationDegreeBlk, "integration_degree_blk"+suffix)
+        gmshFem.userDefinedParameter(_integrationDegreeBlk, "equation_integration_degree_blk"+suffix)
     ))
     {
-      throw Exception("Integration degree (bulk) could not be found.");
+      throw Exception("Equation integration degree (bulk) could not be found.");
     }
     if(!(
-        gmshFem.userDefinedParameter(_integrationDegreeBnd, "integration_degree_bnd")
+        gmshFem.userDefinedParameter(_integrationDegreeBnd, "equation_integration_degree_bnd")
         ||
-        gmshFem.userDefinedParameter(_integrationDegreeBnd, "integration_degree_bnd"+suffix)
+        gmshFem.userDefinedParameter(_integrationDegreeBnd, "equation_integration_degree_bnd"+suffix)
     ))
     {
-      throw Exception("Integration degree (boundary) could not be found.");
+      throw Exception("Equation integration degree (boundary) could not be found.");
     }
     double buffer=0.;
     if(!gmshFem.userDefinedParameter(buffer, "equation_boundary"))
@@ -275,8 +276,9 @@ bool ModelStateConverter::available(Type type) const
 * ParametrizationEquation
 */
 template<Physic T_Physic>
-ParametrizedEquation<T_Physic>::ParametrizedEquation(const ConfigurationInterface* const config, EquationInterface<T_Physic>* const equation, const ParametrizationInterface* const parametrization, const gmshfem::common::GmshFem& gmshFem, std::string suffix) : EquationInterface<T_Physic>(config, gmshFem,suffix), _msc(config), _equation(equation),_parametrization(parametrization)
+ParametrizedEquation<T_Physic>::ParametrizedEquation(const ParametrizationInterface* const parametrization, double pulsation, const ConfigurationInterface* const config, const wave::Discretization<T_Physic>& w_discret, const gmshfem::common::GmshFem& gmshFem, std::string suffix) : EquationInterface<T_Physic>(config,gmshFem,suffix), _msc(config), _equation(newEquation(pulsation,config,w_discret,gmshFem,suffix)),_parametrization(parametrization->clone())
 {
+    _parametrization->pulsation(pulsation);
     if(!_equation->compatible(_parametrization))
     {
         throw Exception("Equation and parametrization are not compatible.");

common/wave/equation/equation.h

@@ -116,10 +116,12 @@ protected:
 private:
     ModelStateConverter _msc;
     EquationInterface<T_Physic>* const _equation;
-    const ParametrizationInterface* const _parametrization;
+    ParametrizationInterface* const _parametrization;
 
 public:
-    ParametrizedEquation(const ConfigurationInterface* const config, EquationInterface<T_Physic>* const equation, const ParametrizationInterface* const parametrization, const gmshfem::common::GmshFem& gmshFem, std::string suffix = "");
+    ParametrizedEquation(const ParametrizationInterface* const parametrization, double pulsation, const ConfigurationInterface* const config, const wave::Discretization<T_Physic>& w_discret, const gmshfem::common::GmshFem& gmshFem, std::string suffix = "");
+
+    ~ParametrizedEquation() {delete _parametrization; delete _equation;};
 
     virtual bool compatible(const ParametrizationInterface* const parametrization) const {return false;};
 

directional.cpp

@@ -12,7 +12,7 @@ using namespace gmshfem::function;
 #include "specific/physic.h"
 #include "specific/configuration/newConfiguration.h"
 #include "specific/model/parametrization/newParametrization.h"
-#include "specific/wave/equation/newEquation.h"
+#include "common/wave/equation/equation.h"
 #include "specific/data/objective/newObjective.h"
 #include "common/statefunctional.h"
 
@@ -46,30 +46,52 @@ SpatialDistribution to_spatialdistribution(const GmshFem& gmshFem)
 template <Physic T_Physic>
 int directional(const GmshFem& gmshFem)
 {
-    std::string name = "noname";
-    gmshFem.userDefinedParameter(name, "name");
+    std::string name = "";
+    gmshFem.userDefinedParameter(name, "prename");
+    std::string suffix = "noname";
+    gmshFem.userDefinedParameter(suffix, "name");
+    name += suffix;
 
     ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
 
     ConfigurationInterface* configuration = newConfiguration(name, parametrization, gmshFem);
 
-    wave::Discretization<T_Physic> w_discret(gmshFem);
-
-    double frequency = 0.;
-    if(!gmshFem.userDefinedParameter(frequency, "frequency0"))
+    unsigned int n_freq = 0;
+    if(!gmshFem.userDefinedParameter(n_freq, "n_freq0"))
     {
-        throw Exception("Frequency could not be found.");
+        throw Exception("Total frequency number could not be found.");
     }
+    std::vector<double> frequency(n_freq);
+    std::vector<const Data<T_Physic>*> d0(n_freq);
+    std::vector<EquationInterface<T_Physic>*> pequation(n_freq);
+    std::vector<ObjectiveInterface<T_Physic>*> objective(n_freq);
 
-    parametrization->pulsation(2.*M_PI*frequency);
-    EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem);
-    EquationInterface<T_Physic>* const pequation = new ParametrizedEquation<T_Physic>(configuration,equation,parametrization,gmshFem);
+    for (unsigned int f = 0; f < n_freq; f++)
+    {
+        std::string suffix_f = std::to_string(f);
+        unsigned int index = 0;
+        if(!gmshFem.userDefinedParameter(index, "frequency0"+suffix_f))
+        {
+            throw Exception("Frequency index #"+suffix_f+" could not be found.");
+        }
+        std::string suffix_idx = std::to_string(index);
 
+        if(!gmshFem.userDefinedParameter(frequency[f], "frequency"+suffix_idx))
+        {
+            throw Exception("Frequency #"+suffix_idx+" could not be found.");
+        }
         std::string filename = "noname";
-    gmshFem.userDefinedParameter(filename, "data0");
-    Data<T_Physic> d0(filename,configuration);
+        if(!gmshFem.userDefinedParameter(filename, "data"+suffix_idx))
+        {
+            throw Exception("Data filename #"+suffix_idx+" could not be found.");
+        }
+        d0[f] = new Data<T_Physic>(filename,configuration);
 
-    ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(&d0,configuration,gmshFem);
+        wave::Discretization<T_Physic> w_discret(gmshFem,suffix_idx);
+        pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,2.*M_PI*frequency[f],configuration,w_discret,gmshFem,suffix_idx);
+
+        objective[f] = newObjective<T_Physic>(d0[f],configuration,gmshFem,suffix_idx);
+    }
 
     StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,nullptr,pequation,objective);
     FunctionalInterface* const functional = statefunctional;
@@ -148,6 +170,8 @@ int directional(const GmshFem& gmshFem)
     msg::print << "Compute directional derivatives" << msg::endl;
     msg::indent();
     Interval interval = to_interval(gmshFem);
+    unsigned int write_data_fields = 0;
+    gmshFem.userDefinedParameter(write_data_fields, "write_data_fields");
     for (unsigned int n = 0; n <= N; n++)
     {
         msg::print << "--- Validation point #" << n << " --- " << msg::endl;
@@ -191,19 +215,23 @@ int directional(const GmshFem& gmshFem)
 
         double j =  functional->performance(m);
         msg::print << "j: "<< j << msg::endl;
-
-        double dj_s = statefunctional->directional(Order::FST, Support::BND, dm);
-        double dj_v = statefunctional->directional(Order::FST, Support::BLK, dm);
+        if(write_data_fields){statefunctional->write_data(Type::F,name+"_data_field"+std::to_string(n));}
 
         /*
+        double dj_s = 0.;
+        double dj_v = 0.;
         double jfd = 0.;
         double d2j_s = 0.;
         double d2j_v = 0.;
         double j2fd = 0.;
         */
 
+        double dj_s = statefunctional->directional(Order::FST, Support::BND, dm);
+        double dj_v = statefunctional->directional(Order::FST, Support::BLK, dm);
         double d2j_s = statefunctional->directional(Order::SCD, Support::BND, dm);
         double d2j_v = statefunctional->directional(Order::SCD, Support::BLK, dm);
+        /*
+        */
 
         double jp =  functional->performance(m + dm*eps);
         //msg::print << "jp: "<< jp << msg::endl;
@@ -225,9 +253,9 @@ int directional(const GmshFem& gmshFem)
         msg::print << "d2j_v: "<< d2j_v << msg::endl;
         msg::print << "d2j: "<< d2j_s+d2j_v << msg::endl;
         msg::print << "j2FD: "<< j2fd << msg::endl;
-
         /*
         */
+
         for (unsigned int c = 0; c < parametrization->size(); c++)
         {
             output << mn[c];
@@ -241,9 +269,11 @@ int directional(const GmshFem& gmshFem)
     output.close();
 
     delete statefunctional;
-    delete equation;
-    delete pequation;
-    delete objective;
+    for (unsigned int f = 0; f < n_freq; f++)
+    {
+        delete pequation[f];
+        delete objective[f];
+    }
     delete configuration;
     delete parametrization;
 

gradient.cpp

@@ -12,7 +12,7 @@ using namespace gmshfem::common;
 #include "specific/physic.h"
 #include "specific/configuration/newConfiguration.h"
 #include "specific/model/parametrization/newParametrization.h"
-#include "specific/wave/equation/newEquation.h"
+#include "common/wave/equation/equation.h"
 #include "specific/data/objective/newObjective.h"
 #include "specific/model/innerproduct/sobolev.h"
 #include "common/statefunctional.h"
@@ -40,35 +40,60 @@ int gradient(const GmshFem& gmshFem)
     gmshFem.userDefinedParameter(suffix, "name");
     name += suffix;
 
-    ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
+    const ParametrizationInterface* const parametrization = newParametrization<T_Physic>(gmshFem);
 
-    ConfigurationInterface* configuration = newConfiguration(name, parametrization, gmshFem);
+    const ConfigurationInterface* const configuration = newConfiguration(name, parametrization, gmshFem);
 
 
-    wave::Discretization<T_Physic> w_discret(gmshFem,"0");
-    model::Discretization m_discret(gmshFem,"0");
+    model::Discretization m_discret(gmshFem);
 
-    double frequency = 0.;
-    if(!gmshFem.userDefinedParameter(frequency, "frequency0"))
+    unsigned int n_freq = 0;
+    if(!gmshFem.userDefinedParameter(n_freq, "n_freq0"))
     {
-        throw Exception("Frequency could not be found.");
+        throw Exception("Total frequency number could not be found.");
     }
+    std::vector<double> frequency(n_freq);
+    std::vector<const Data<T_Physic>*> d0(n_freq);
+    std::vector<EquationInterface<T_Physic>*> pequation(n_freq);
+    std::vector<ObjectiveInterface<T_Physic>*> objective(n_freq);
 
-    parametrization->pulsation(2.*M_PI*frequency);
-    EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem);
-    EquationInterface<T_Physic>* const pequation = new ParametrizedEquation<T_Physic>(configuration,equation,parametrization,gmshFem);
+    for (unsigned int f = 0; f < n_freq; f++)
+    {
+        std::string suffix_f = std::to_string(f);
+        unsigned int index = 0;
+        if(!gmshFem.userDefinedParameter(index, "frequency0"+suffix_f))
+        {
+            throw Exception("Frequency index #"+suffix_f+" could not be found.");
+        }
+        std::string suffix_idx = std::to_string(index);
 
+        if(!gmshFem.userDefinedParameter(frequency[f], "frequency"+suffix_idx))
+        {
+            throw Exception("Frequency #"+suffix_idx+" could not be found.");
+        }
         std::string filename = "noname";
-    gmshFem.userDefinedParameter(filename, "data0");
-    Data<T_Physic> d0(filename,configuration);
+        if(!gmshFem.userDefinedParameter(filename, "data"+suffix_idx))
+        {
+            throw Exception("Data filename #"+suffix_idx+" could not be found.");
+        }
+        d0[f] = new Data<T_Physic>(filename,configuration);
 
-    ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(&d0,configuration,gmshFem);
+        wave::Discretization<T_Physic> w_discret(gmshFem,suffix_idx);
+        pequation[f] = new ParametrizedEquation<T_Physic>(parametrization,2.*M_PI*frequency[f],configuration,w_discret,gmshFem,suffix_idx);
 
-    sobolev::InnerProduct innerproduct(configuration,m_discret,gmshFem);
+        objective[f] = newObjective<T_Physic>(d0[f],configuration,gmshFem,suffix_idx);
+    }
+
+    sobolev_l2::InnerProduct innerproduct(configuration,m_discret,gmshFem);
 
     StateFunctional<T_Physic>* const statefunctional = new StateFunctional<T_Physic>(configuration,&innerproduct,pequation,objective);
     FunctionalInterface* const functional = statefunctional;
-    functional->setModel(configuration->m0(2.*M_PI*frequency));
+    functional->setModel(configuration->m0(2.*M_PI*frequency[0]));
+    for (unsigned int c = 0; c < parametrization->size(); c++)
+    {
+        gmshfem::post::save(configuration->m0(2.*M_PI*frequency[0])[c], configuration->model_unknown(Support::BLK), name+"_m_blkc"+std::to_string(c)+"f0", "pos", "");
+        gmshfem::post::save(configuration->m0(2.*M_PI*frequency[0])[c], configuration->model_unknown(Support::BND), name+"_m_bndc"+std::to_string(c)+"f0", "pos", "");
+    }
 
     double scale0,scaleN;
     unsigned int N;
@@ -106,7 +131,7 @@ int gradient(const GmshFem& gmshFem)
         msg::print << "--- Scale #"<< n << " = "<< scale << " --- " << msg::endl;
         msg::indent();
 
-        innerproduct.alpha1( (scale / (2 * M_PI)) * (scale / (2 * M_PI)) );
+        innerproduct.alpha( (scale / (2 * M_PI)) * (scale / (2 * M_PI)) );
         innerproduct.bulkWeight(1.0);
         innerproduct.boundaryWeight(0.0);
         statefunctional->innerproductIsObsolete();
@@ -122,9 +147,11 @@ int gradient(const GmshFem& gmshFem)
     msg::unindent();
 
     delete statefunctional;
-    delete equation;
-    delete pequation;
-    delete objective;
+    for (unsigned int f = 0; f < n_freq; f++)
+    {
+        delete pequation[f];
+        delete objective[f];
+    }
     delete configuration;
     delete parametrization;
 

input/Bones1D/common.txt

+prename = bones1d_
+#Configuration
+configuration = line_acquisition
+#
+#Parametrization
+parametrization = natural
+model_reference_frequency = 1.
+#
+receiver_on_emitter = 1
+xer = 0.5
+L = 6.
+#
+Re(mbc0) = 1.
+Im(mbc0) = 0.
+#
+#Equation
+physic = acoustic
+equation = helmholtz
+equation_boundary = 1
+#
+#Discretization
+h = 0.01
+integration_type = Gauss
+#Wave
+wave_FunctionSpaceType = HierarchicalH1
+wave_FunctionSpaceDegree = 1
+#
+equation_integration_degree_bnd = 0
+equation_integration_degree_blk = 3
+innerproduct_integration_degree_blk = 2
+innerproduct_integration_degree_bnd = 0
+#
+n_freq=20
+#Frequency
+frequency0 = 0.125
+frequency1 = 0.25
+frequency2 = 0.375
+frequency3 = 0.5
+frequency4 = 0.625
+frequency5 = 0.75
+frequency6 = 0.875
+frequency7 = 1.0
+frequency8 = 1.25
+frequency9 = 1.5
+frequency10 = 1.75
+frequency11 = 2.0
+frequency12 = 2.25
+frequency13 = 2.5
+frequency14 = 2.75
+frequency15 = 3.0
+frequency16 = 3.25
+frequency17 = 3.5
+frequency18 = 3.75
+frequency19 = 4.0

input/Bones1D/synthetic.txt

+name = synthetics
+#Configuration
+unknown = none
+inclusion_ni=3
+#Inclusion 0
+inclusion_position0 = 2.
+inclusion_depth0 = 1.0
+Re(mi0c0) = 1.2
+Im(mi0c0) = 0.
+#Inclusion 1
+inclusion_position1 = 3.0
+inclusion_depth1 = 1.0
+Re(mi1c0) = 0.8
+Im(mi1c0) = 0.
+#Inclusion 2
+inclusion_position2 = 4.0
+inclusion_depth2 = 1.0
+Re(mi2c0) = 1.3
+Im(mi2c0) = 0.
+#
+#Output
+write_model_fields = 0
+write_wave_fields = 0
+write_data_fields = 0

input/EM2007/ModelX/cluster_syntetics.sh

@@ -4,7 +4,7 @@
 #SBATCH --output=EM2007ModelX_synthetics.log
 
 #SBATCH --ntasks=1
-#SBATCH --time=30:00
+#SBATCH --time=60:00
 #SBATCH --mem-per-cpu=48000
 
 time ./synthetics ../input/EM2007/ModelX/common.txt ../input/EM2007/ModelX/synthetics.txt -verbose 2 -maxThreads 1

input/EM2007/ModelX/common.txt

@@ -57,35 +57,35 @@ n_freq = 7
 #Frequency0
 frequency0 = 0.25
 wave_FunctionSpaceDegree0 = 1
-integration_degree_blk0 = 3
-integration_degree_bnd0 = 4
+equation_integration_degree_blk0 = 3
+equation_integration_degree_bnd0 = 4
 #Frequency1
 frequency1 = 0.3
 wave_FunctionSpaceDegree1 = 2
-integration_degree_blk1 = 5
-integration_degree_bnd1 = 6
+equation_integration_degree_blk1 = 5
+equation_integration_degree_bnd1 = 6
 #Frequency2
 frequency2 = 0.4
 wave_FunctionSpaceDegree2 = 2
-integration_degree_blk2 = 5
-integration_degree_bnd2 = 6
+equation_integration_degree_blk2 = 5
+equation_integration_degree_bnd2 = 6
 #Frequency3
 frequency3 = 0.5
 wave_FunctionSpaceDegree3 = 2
-integration_degree_blk3 = 5
-integration_degree_bnd3 = 6
+equation_integration_degree_blk3 = 5
+equation_integration_degree_bnd3 = 6
 #Frequency4
 frequency4 = 0.6
 wave_FunctionSpaceDegree4 = 3
-integration_degree_blk4 = 7
-integration_degree_bnd4 = 8
+equation_integration_degree_blk4 = 7
+equation_integration_degree_bnd4 = 8
 #Frequency5
 frequency5 = 0.75
 wave_FunctionSpaceDegree5 = 4
-integration_degree_blk5 = 9
-integration_degree_bnd5 = 10
+equation_integration_degree_blk5 = 9
+equation_integration_degree_bnd5 = 10
 #Frequency6
 frequency6 = 1.0
 wave_FunctionSpaceDegree6 = 5
-integration_degree_blk6 = 11
-integration_degree_bnd6 = 12
+equation_integration_degree_blk6 = 11
+equation_integration_degree_bnd6 = 12

input/EM2007/ModelX/inversion_nscgB.txt

 #
 #Local minimum search
 localminimumsearch = yuanfan
-localminimum_initMu = 1
+localminimum_retrospective = 0
+localminimum_initMu = 0.25
 localminimum_c0 = 0.0001
 localminimum_c1 = 0.75
 localminimum_c2 = 0.75
 localminimum_c3 = 0.25
 localminimum_c4 = 2.
-localminimum_writeInterval = 5
-localminimum_maxIteration = 50
+localminimum_writeInterval = 2
+localminimum_maxIteration = 20
 #
 #Descent search
 descentsearch = newton_steihaug_conjugategradient

input/EM2007/ModelX/simultaneous/cluster_inversion.sh

+#!/bin/bash
+#
+#SBATCH --job-name=EM2007ModelX_inversion_sim
+#SBATCH --output=EM2007ModelX_inversion_sim.log
+
+#SBATCH --ntasks=1
+#SBATCH --time=800:00
+#SBATCH --mem-per-cpu=48000
+
+time ./inversion ../input/EM2007/ModelX/common.txt ../input/EM2007/ModelX/simultaneous/inversion.txt ../input/EM2007/ModelX/inversion_nscgB.txt -verbose 2 -maxThreads 1