From c8bd743d44c88acb4c9b9bc9c930227d39ad06d8 Mon Sep 17 00:00:00 2001 From: XavierAdriaens <xavier.adriaens@doct.uliege.be> Date: Wed, 12 Aug 2020 10:24:37 +0200 Subject: [PATCH] reproducing Figure 1 of Virieux & Operto 2009 review paper --- CONVENTIONS.txt | 4 +- common/configuration.h | 5 +- common/functional.h | 3 + common/model/discretization.cpp | 19 ++- common/model/discretization.h | 2 +- common/model/element.cpp | 19 +++ common/model/element.h | 3 + common/model/innerproduct/innerproduct.cpp | 10 +- common/model/innerproduct/innerproduct.h | 6 +- common/optimization/descentsearch.h | 3 +- common/optimization/linesearch.h | 2 +- common/optimization/minimumsearch.h | 2 +- common/rediscretize.cpp | 23 +++ common/rediscretize.h | 10 ++ common/statefunctional.h | 2 + convergence.cpp | 8 +- directional.cpp | 6 +- gradient.cpp | 6 +- input/VO2009Algorithm/common.txt | 46 ++++++ input/VO2009Algorithm/inversion_gd.txt | 39 +++++ input/VO2009Algorithm/inversion_ncg.txt | 41 ++++++ input/VO2009Algorithm/synthetics.txt | 13 ++ input/paper1/common.txt | 2 +- .../paper1/convergence_highly_background.txt | 2 +- .../paper1/convergence_weakly_background.txt | 2 +- .../paper1/directional_highly_background.txt | 2 +- input/paper1/directional_highly_inclusion.txt | 2 +- .../paper1/directional_weakly_background.txt | 2 +- input/paper1/directional_weakly_inclusion.txt | 2 +- input/paper1/gradient_highly.txt | 2 +- inversion.cpp | 112 ++++++++++++++ specific/configuration/inclusion.cpp | 19 ++- specific/configuration/inclusion.h | 7 +- specific/configuration/newConfiguration.h | 6 +- specific/configuration/soil.cpp | 4 +- specific/configuration/soil.h | 3 +- specific/model/innerproduct/newInnerProduct.h | 4 +- specific/model/innerproduct/sobolev.cpp | 19 ++- specific/model/innerproduct/sobolev.h | 2 +- .../descentsearch/gradientdescent.cpp | 18 +++ .../descentsearch/gradientdescent.h | 40 +++++ .../descentsearch/newDescentSearch.h | 31 ++++ .../newton_conjugategradient.cpp | 139 ++++++++++++++++++ .../descentsearch/newton_conjugategradient.h | 72 +++++++++ .../optimization/linesearch/backtracking.cpp | 100 +++++++++++++ .../optimization/linesearch/backtracking.h | 57 +++++++ specific/optimization/linesearch/constant.cpp | 73 +++++++++ specific/optimization/linesearch/constant.h | 53 +++++++ .../optimization/linesearch/newLineSearch.h | 33 +++++ .../optimization/linesearch/secondorder.cpp | 72 +++++++++ .../optimization/linesearch/secondorder.h | 53 +++++++ .../optimization/minimumsearch/classic.cpp | 13 +- specific/optimization/minimumsearch/classic.h | 2 +- specific/wave/discretization.cpp | 19 ++- specific/wave/discretization.h | 2 +- specific/wave/equation/helmholtz.cpp | 54 +++++-- specific/wave/equation/helmholtz.h | 8 +- synthetics.cpp | 39 +++-- 58 files changed, 1243 insertions(+), 99 deletions(-) create mode 100644 common/rediscretize.cpp create mode 100644 common/rediscretize.h create mode 100644 input/VO2009Algorithm/common.txt create mode 100644 input/VO2009Algorithm/inversion_gd.txt create mode 100644 input/VO2009Algorithm/inversion_ncg.txt create mode 100644 input/VO2009Algorithm/synthetics.txt create mode 100644 inversion.cpp create mode 100644 specific/optimization/descentsearch/gradientdescent.cpp create mode 100644 specific/optimization/descentsearch/gradientdescent.h create mode 100644 specific/optimization/descentsearch/newDescentSearch.h create mode 100644 specific/optimization/descentsearch/newton_conjugategradient.cpp create mode 100644 specific/optimization/descentsearch/newton_conjugategradient.h create mode 100644 specific/optimization/linesearch/backtracking.cpp create mode 100644 specific/optimization/linesearch/backtracking.h create mode 100644 specific/optimization/linesearch/constant.cpp create mode 100644 specific/optimization/linesearch/constant.h create mode 100644 specific/optimization/linesearch/newLineSearch.h create mode 100644 specific/optimization/linesearch/secondorder.cpp create mode 100644 specific/optimization/linesearch/secondorder.h diff --git a/CONVENTIONS.txt b/CONVENTIONS.txt index 087edd7..c3928b3 100644 --- a/CONVENTIONS.txt +++ b/CONVENTIONS.txt @@ -1,10 +1,10 @@ /* * 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. /* * Implementation */ -*Passing through references is used only to avoid copy. Passing through pointers is used otherwise. +*Passing through references is used only to avoid copy. Passing through pointers is used otherwise. diff --git a/common/configuration.h b/common/configuration.h index aa834bc..a88c7b6 100644 --- a/common/configuration.h +++ b/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.;}; diff --git a/common/functional.h b/common/functional.h index 8009831..e0da9cb 100644 --- a/common/functional.h +++ b/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; diff --git a/common/model/discretization.cpp b/common/model/discretization.cpp index 6ad71b5..a29f7c9 100644 --- a/common/model/discretization.cpp +++ b/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."); } diff --git a/common/model/discretization.h b/common/model/discretization.h index 463df91..697dcfc 100644 --- a/common/model/discretization.h +++ b/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;}; diff --git a/common/model/element.cpp b/common/model/element.cpp index b57fdef..49e02cd 100644 --- a/common/model/element.cpp +++ b/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 */ diff --git a/common/model/element.h b/common/model/element.h index 9c8f813..00a6218 100644 --- a/common/model/element.h +++ b/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(); }; /* diff --git a/common/model/innerproduct/innerproduct.cpp b/common/model/innerproduct/innerproduct.cpp index dec1536..b7bc12f 100644 --- a/common/model/innerproduct/innerproduct.cpp +++ b/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; diff --git a/common/model/innerproduct/innerproduct.h b/common/model/innerproduct/innerproduct.h index 92cb5b7..4f77a56 100644 --- a/common/model/innerproduct/innerproduct.h +++ b/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; diff --git a/common/optimization/descentsearch.h b/common/optimization/descentsearch.h index 64de571..86f8a62 100644 --- a/common/optimization/descentsearch.h +++ b/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; }; diff --git a/common/optimization/linesearch.h b/common/optimization/linesearch.h index b9ad1d9..df2edfd 100644 --- a/common/optimization/linesearch.h +++ b/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; }; diff --git a/common/optimization/minimumsearch.h b/common/optimization/minimumsearch.h index e3bccc1..d4ac180 100644 --- a/common/optimization/minimumsearch.h +++ b/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; }; diff --git a/common/rediscretize.cpp b/common/rediscretize.cpp new file mode 100644 index 0000000..40f064f --- /dev/null +++ b/common/rediscretize.cpp @@ -0,0 +1,23 @@ +//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; +} diff --git a/common/rediscretize.h b/common/rediscretize.h new file mode 100644 index 0000000..ed7d440 --- /dev/null +++ b/common/rediscretize.h @@ -0,0 +1,10 @@ +#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 diff --git a/common/statefunctional.h b/common/statefunctional.h index 7834c8f..a6fe184 100644 --- a/common/statefunctional.h +++ b/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); diff --git a/convergence.cpp b/convergence.cpp index 2f6d816..6246040 100644 --- a/convergence.cpp +++ b/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); @@ -83,7 +83,7 @@ int convergence(const GmshFem& gmshFem) double j = functional->performance( m ); msg::print << "j: "<< j << msg::endl; - + double dj_v = statefunctional->directional(Order::FST, Support::BLK, dm); msg::print << "dj_v: "<< dj_v << msg::endl; diff --git a/directional.cpp b/directional.cpp index 9364d80..4c93512 100644 --- a/directional.cpp +++ b/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); diff --git a/gradient.cpp b/gradient.cpp index 9c3586f..fdb8434 100644 --- a/gradient.cpp +++ b/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); diff --git a/input/VO2009Algorithm/common.txt b/input/VO2009Algorithm/common.txt new file mode 100644 index 0000000..6660fc0 --- /dev/null +++ b/input/VO2009Algorithm/common.txt @@ -0,0 +1,46 @@ +#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 diff --git a/input/VO2009Algorithm/inversion_gd.txt b/input/VO2009Algorithm/inversion_gd.txt new file mode 100644 index 0000000..0a18cfc --- /dev/null +++ b/input/VO2009Algorithm/inversion_gd.txt @@ -0,0 +1,39 @@ +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 diff --git a/input/VO2009Algorithm/inversion_ncg.txt b/input/VO2009Algorithm/inversion_ncg.txt new file mode 100644 index 0000000..a5e27d8 --- /dev/null +++ b/input/VO2009Algorithm/inversion_ncg.txt @@ -0,0 +1,41 @@ +name = VO2009A_ncg +#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 = 0 +# +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 = newton_conjugategradient +descent_maxIteration = 3 +eta20 = 0.0 +# +#Line search +linesearch = constant +a = 1 diff --git a/input/VO2009Algorithm/synthetics.txt b/input/VO2009Algorithm/synthetics.txt new file mode 100644 index 0000000..ab5ffa9 --- /dev/null +++ b/input/VO2009Algorithm/synthetics.txt @@ -0,0 +1,13 @@ +name = VO2009A +#Configuration +inclusion = cylinder +filled = 1 +unknown = none +# +r = 0.2 +# +Re(mc) = 0.9 +Im(mc) = 0. +# +#Discretization +h = 0.05 diff --git a/input/paper1/common.txt b/input/paper1/common.txt index abe244a..1d801e5 100644 --- a/input/paper1/common.txt +++ b/input/paper1/common.txt @@ -25,4 +25,4 @@ wave_IntegrationType = Gauss wave_FunctionSpaceDegree = 5 # #Frequency -frequency = 1. +frequency0 = 1. diff --git a/input/paper1/convergence_highly_background.txt b/input/paper1/convergence_highly_background.txt index 4e1e376..da90e33 100644 --- a/input/paper1/convergence_highly_background.txt +++ b/input/paper1/convergence_highly_background.txt @@ -23,7 +23,7 @@ epsN = 1e-1 N = 9 # #Data -data = filled_inclusion_data +data0 = filled_inclusion_data0 # #Discretization h = 0.25 diff --git a/input/paper1/convergence_weakly_background.txt b/input/paper1/convergence_weakly_background.txt index 1823a5a..0d59cbf 100644 --- a/input/paper1/convergence_weakly_background.txt +++ b/input/paper1/convergence_weakly_background.txt @@ -23,7 +23,7 @@ epsN = 1e-1 N = 9 # #Data -data = empty_inclusion_data +data0 = empty_inclusion_data0 # #Discretization h = 0.25 diff --git a/input/paper1/directional_highly_background.txt b/input/paper1/directional_highly_background.txt index bf08322..2af2d10 100644 --- a/input/paper1/directional_highly_background.txt +++ b/input/paper1/directional_highly_background.txt @@ -25,7 +25,7 @@ N = 20 eps = 1e-5 # #Data -data = filled_inclusion_data +data0 = filled_inclusion_data0 # #Discretization h = 0.25 diff --git a/input/paper1/directional_highly_inclusion.txt b/input/paper1/directional_highly_inclusion.txt index 0883011..c38d407 100644 --- a/input/paper1/directional_highly_inclusion.txt +++ b/input/paper1/directional_highly_inclusion.txt @@ -25,7 +25,7 @@ N = 20 eps = 1e-5 # #Data -data = filled_inclusion_data +data0 = filled_inclusion_data0 # #Discretization h = 0.25 diff --git a/input/paper1/directional_weakly_background.txt b/input/paper1/directional_weakly_background.txt index 1a90153..6a29097 100644 --- a/input/paper1/directional_weakly_background.txt +++ b/input/paper1/directional_weakly_background.txt @@ -25,7 +25,7 @@ N = 20 eps = 1e-5 # #Data -data = empty_inclusion_data +data0 = empty_inclusion_data0 # #Discretization h = 0.25 diff --git a/input/paper1/directional_weakly_inclusion.txt b/input/paper1/directional_weakly_inclusion.txt index ae23fbd..ba98a57 100644 --- a/input/paper1/directional_weakly_inclusion.txt +++ b/input/paper1/directional_weakly_inclusion.txt @@ -26,7 +26,7 @@ N = 20 eps = 1e-5 # #Data -data = empty_inclusion_data +data0 = empty_inclusion_data0 # #Discretization h = 0.25 diff --git a/input/paper1/gradient_highly.txt b/input/paper1/gradient_highly.txt index c21394e..4bd3a51 100644 --- a/input/paper1/gradient_highly.txt +++ b/input/paper1/gradient_highly.txt @@ -17,7 +17,7 @@ interval = log N = 2 # #Data -data = filled_inclusion_data +data0 = filled_inclusion_data0 # #Inner product complex = 0 diff --git a/inversion.cpp b/inversion.cpp new file mode 100644 index 0000000..1e19659 --- /dev/null +++ b/inversion.cpp @@ -0,0 +1,112 @@ +//GmshFem Library +#include "GmshFem.h" + +using namespace gmshfem; +using namespace gmshfem::common; + +//FWI Library +#include "specific/physic.h" +#include "specific/configuration/newConfiguration.h" +#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/statefunctional.h" +#include "specific/optimization/minimumsearch/newMinimumSearch.h" +#include "specific/optimization/descentsearch/newDescentSearch.h" +#include "specific/optimization/linesearch/newLineSearch.h" + +template <Physic T_Physic> +int inversion(const GmshFem& gmshFem) +{ + std::string name = "noname"; + gmshFem.userDefinedParameter(name, "name"); + + ConfigurationInterface* configuration = newConfiguration(name, 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); + + unsigned int n_group = 0; + if(!gmshFem.userDefinedParameter(n_group, "n_group")) + { + throw Exception("Number of frequency group could not be found."); + } + + ModelField m; + for (unsigned int g = 0; g < n_group; g++) + { + std::string suffix = std::to_string(g); + msg::print << "Inverse frequency group #" + suffix << msg::endl; + std::string integrationType; + gmshFem.userDefinedParameter(integrationType, "model_IntegrationType"); + gmshFem.userDefinedParameter(integrationType, "model_IntegrationType"+suffix); + model::Discretization m_discret(gmshFem,suffix); + if(g==0) + { + m = initialize(configuration, m_discret, integrationType); + } + else + { + m = rediscretize(m, configuration, m_discret, integrationType); + } + m.write(name+"_m"+suffix); + + double frequency; + if(!gmshFem.userDefinedParameter(frequency, "frequency"+suffix)) + { + throw Exception("Frequency "+suffix+" could not be found."); + } + + wave::Discretization<T_Physic> w_discret(gmshFem,suffix); + EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem); + + std::string filename = "noname"; + gmshFem.userDefinedParameter(filename, "data"+suffix); + Data<T_Physic> d0(filename,configuration); + + ObjectiveInterface<T_Physic>* const objective = newObjective<T_Physic>(&d0,gmshFem); + + InnerProductInterface* const innerproduct = newInnerProduct(complex,configuration,m_discret,gmshFem,suffix); + + FunctionalInterface* const functional = new StateFunctional<T_Physic>(configuration,innerproduct,equation,objective); + + MinimumSearchInterface* const minimumsearch = newMinimumSearch(gmshFem); + DescentSearchInterface* const descentsearch = newDescentSearch(gmshFem); + LineSearchInterface* const linesearch = newLineSearch(gmshFem); + + (*minimumsearch)(&m, functional, descentsearch, linesearch); + minimumsearch->history()->write(name+"_minimum_history" + suffix); + descentsearch->history()->write(name+"_descent_history" + suffix); + linesearch->history()->write(name+"_linesearch_history" + suffix); + + delete minimumsearch; + delete descentsearch; + delete linesearch; + delete functional; + delete innerproduct; + delete objective; + delete equation; + } + std::string suffix = std::to_string(n_group); + m.write(name+"_m"+suffix); + + delete configuration; + return 0; +} + +int main(int argc, char **argv) +{ + GmshFem gmshFem(argc, argv); + Physic T_Physic = to_physic(gmshFem); + switch (T_Physic) + { + case Physic::acoustic: default: + return inversion<Physic::acoustic>(gmshFem); + } +} diff --git a/specific/configuration/inclusion.cpp b/specific/configuration/inclusion.cpp index 7bf3891..18b3040 100644 --- a/specific/configuration/inclusion.cpp +++ b/specific/configuration/inclusion.cpp @@ -56,7 +56,7 @@ namespace inclusion /* * class Configuration */ - Configuration::Configuration(const GmshFem& gmshFem) : ConfigurationInterface(gmshFem) + Configuration::Configuration(std::string name, const GmshFem& gmshFem) : ConfigurationInterface(name, gmshFem) { msg::print << "Initialize inclusion configuration." << msg::endl; @@ -253,7 +253,7 @@ namespace inclusion (_inclusion_type!=Inclusion::None) ) { - throw common::Exception("Inclusion model parameter could not be found."); + throw Exception("Inclusion model parameter could not be found."); } _mc = Remc + im * Immc; @@ -268,7 +268,7 @@ namespace inclusion msg::print << "Generate meshes" << msg::endl; gmsh::option::setNumber("General.Terminal", 1); - std::string name = "inclusion"; + std::string name = _name + "_inclusion"; switch (_inclusion_type) { case Inclusion::Cylinder: name += "_cyl"; break; @@ -448,6 +448,7 @@ namespace inclusion gmodel::setPhysicalName(1, 13, "inclusion_bnd"); } } + factory::removeAllDuplicates(); } void Configuration::data_mesh() const { @@ -511,7 +512,17 @@ namespace inclusion template<Physic T_Physic> ModelFunction green0_preconditioner(const WaveMultiField<T_Physic>& g, const Configuration* const config) { - return 1.; + if( !config->receiversAreEmittersX() && !config->receiversAreEmittersY() ) + { + return autocorrelate<T_Physic>(g,config->emitter_idx_X())*autocorrelate<T_Physic>(g,config->receiver_idx_X()) + + + autocorrelate<T_Physic>(g,config->emitter_idx_Y())*autocorrelate<T_Physic>(g,config->receiver_idx_Y()); + } + else + { + throw Exception("Model preconditioner not implemented."); + return 0.; + } } template ModelFunction green0_preconditioner<Physic::acoustic>(const WaveMultiField<Physic::acoustic>&, const Configuration* const); diff --git a/specific/configuration/inclusion.h b/specific/configuration/inclusion.h index 1d7a391..ec4a0e0 100644 --- a/specific/configuration/inclusion.h +++ b/specific/configuration/inclusion.h @@ -69,8 +69,12 @@ namespace inclusion virtual void wave_mesh() const; virtual void data_mesh() const; public: - Configuration(const gmshfem::common::GmshFem& gmshFem); + Configuration(std::string name, const gmshfem::common::GmshFem& gmshFem); + const std::vector<unsigned int>& emitter_idx_X() const {return _emitter_idx_X;}; + const std::vector<unsigned int>& emitter_idx_Y() const {return _emitter_idx_Y;}; + const std::vector<unsigned int>& receiver_idx_X() const {return _receiver_idx_X;}; + const std::vector<unsigned int>& receiver_idx_Y() const {return _receiver_idx_Y;}; bool receiversAreEmittersX () const {return (std::abs(_xe-_xr)<1e-14) && (std::abs(_He-_Hr)<1e-14) && (_nxe==_nxr);}; bool receiversAreEmittersY () const {return (std::abs(_ye-_yr)<1e-14) && (std::abs(_Le-_Lr)<1e-14) && (_nye==_nyr);}; bool isFilled() const {return _isFilled;}; @@ -81,6 +85,7 @@ namespace inclusion unsigned int nyr() const {return _nyr;}; virtual void mesh() const; + virtual double area() const {return _L * _H; }; 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);}; diff --git a/specific/configuration/newConfiguration.h b/specific/configuration/newConfiguration.h index 3797abe..25eaedc 100644 --- a/specific/configuration/newConfiguration.h +++ b/specific/configuration/newConfiguration.h @@ -12,15 +12,15 @@ #include "soil.h" #include "inclusion.h" -ConfigurationInterface* newConfiguration(const gmshfem::common::GmshFem& gmshFem) +ConfigurationInterface* newConfiguration(std::string name, const gmshfem::common::GmshFem& gmshFem) { std::string configuration; if(!gmshFem.userDefinedParameter(configuration, "configuration")) { throw gmshfem::common::Exception("Configuration type could not be found."); } - if(configuration=="soil"){return new soil::Configuration(gmshFem);} - else if(configuration=="inclusion"){return new inclusion::Configuration(gmshFem);} + if(configuration=="soil"){return new soil::Configuration(name, gmshFem);} + else if(configuration=="inclusion"){return new inclusion::Configuration(name, gmshFem);} else { throw gmshfem::common::Exception("Configuration " + configuration + " is not valid."); diff --git a/specific/configuration/soil.cpp b/specific/configuration/soil.cpp index d6754b9..83d4283 100644 --- a/specific/configuration/soil.cpp +++ b/specific/configuration/soil.cpp @@ -110,7 +110,7 @@ namespace soil /* * class Configuration */ - Configuration::Configuration(const GmshFem& gmshFem) : ConfigurationInterface(gmshFem) + Configuration::Configuration(std::string name, const GmshFem& gmshFem) : ConfigurationInterface(name, gmshFem) { msg::print << "Initialize soil configuration." << msg::endl; @@ -250,7 +250,7 @@ namespace soil { msg::print << "Generate meshes" << msg::endl; gmsh::option::setNumber("General.Terminal", 1); - std::string name = "soil"; + std::string name = _name + "_soil"; gmodel::add("soil"); if(!_remesh) { diff --git a/specific/configuration/soil.h b/specific/configuration/soil.h index 5e00d4a..2af27c4 100644 --- a/specific/configuration/soil.h +++ b/specific/configuration/soil.h @@ -41,10 +41,11 @@ namespace soil virtual void wave_mesh() const; virtual void data_mesh() const; public: - Configuration(const gmshfem::common::GmshFem& gmshFem); + Configuration(std::string name, const gmshfem::common::GmshFem& gmshFem); virtual void mesh() const; unsigned int ner() const {return _ner;}; double H() const {return _H;}; + virtual double area() const {return _L * _H;}; virtual double array() const {return _Ler;}; virtual double depth() const {return _H;}; virtual double m_reference() const {return std::real(_m_soil);}; diff --git a/specific/model/innerproduct/newInnerProduct.h b/specific/model/innerproduct/newInnerProduct.h index 49d7a6b..587a28c 100644 --- a/specific/model/innerproduct/newInnerProduct.h +++ b/specific/model/innerproduct/newInnerProduct.h @@ -11,14 +11,14 @@ //GmshFWI Library #include "sobolev.h" -InnerProductInterface* newInnerProduct(bool complex, 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 suffix) { 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(complex,config,m_discret,gmshFem);} + if(innerproduct=="sobolev"){return new sobolev::InnerProduct(complex,config,m_discret,gmshFem, suffix);} else { throw gmshfem::common::Exception("Inner product " + innerproduct + " is not valid."); diff --git a/specific/model/innerproduct/sobolev.cpp b/specific/model/innerproduct/sobolev.cpp index dc3757c..d9f14fd 100644 --- a/specific/model/innerproduct/sobolev.cpp +++ b/specific/model/innerproduct/sobolev.cpp @@ -17,6 +17,13 @@ static const std::complex< double > im = std::complex< double >(0., 1.); */ namespace sobolev { + InnerProduct::InnerProduct(bool complex, const ConfigurationInterface* const config, const model::Discretization& m_discret, const GmshFem& gmshFem, std::string suffix) : DifferentialInnerProductInterface(complex, config,m_discret,gmshFem,suffix), _mu(nullptr), _su(nullptr), _diag_preconditioner(false), _alpha1(0.), _bulkWeight(1.0), _boundaryWeight(0.) + { + unsigned int diag_preconditioner = 0; + gmshFem.userDefinedParameter(diag_preconditioner, "diag_preconditioner"); + _diag_preconditioner = ((bool) diag_preconditioner); + }; + void InnerProduct::link(ModelUpdater* const mu, SensitivityUpdater* const su) { if(_mu!=nullptr || _su!=nullptr) @@ -48,14 +55,16 @@ namespace sobolev _formulation.galerkin(_alpha1*grad(dof(_j)), grad(tf(_j)), _config->model_unknown(Support::BLK),integrationType(2*degree)); std::array<bool,4> needToBeUpToDate = {true,false,false,false}; - + + //gmshfem::post::save(_su->get(Order::DIAG,Support::BLK,_mu->get(needToBeUpToDate)),_config->model_unknown(Support::BLK),"precond","pos",""); + //throw Exception("STOOOOOOOOOOOOOOOOOOP."); if(_diag_preconditioner) { - _formulation.galerkin( _su->get(Order::DIAG,Support::BLK,_mu->get(needToBeUpToDate)) * dof(_j), tf(_j), _config->model_unknown(Support::BLK), integrationType(3*degree)); + _formulation.galerkin(_su->get(Order::DIAG,Support::BLK,_mu->get(needToBeUpToDate)) * dof(_j), tf(_j), _config->model_unknown(Support::BLK), integrationType(3*degree)); } else { - _formulation.galerkin( dof(_j), tf(_j), _config->model_unknown(Support::BLK), integrationType(3*degree)); + _formulation.galerkin(dof(_j), tf(_j), _config->model_unknown(Support::BLK), integrationType(3*degree)); } } void InnerProduct::setRHS(const ModelFunction& bulk_sensi,const ModelFunction& boundary_sensi) @@ -64,11 +73,11 @@ namespace sobolev if(_bulkWeight!=0.) { - _formulation.galerkin(-_bulkWeight * bulk_sensi,tf(_j), _config->model_unknown(Support::BLK), integrationType(2*degree)); + _formulation.galerkin(-_bulkWeight * conj(bulk_sensi),tf(_j), _config->model_unknown(Support::BLK), integrationType(2*degree)); } if(_boundaryWeight!=0.) { - _formulation.galerkin(-_boundaryWeight * boundary_sensi,tf(_j), _config->model_unknown(Support::BND), integrationType(2*degree)); + _formulation.galerkin(-_boundaryWeight * conj(boundary_sensi),tf(_j), _config->model_unknown(Support::BND), integrationType(2*degree)); } } diff --git a/specific/model/innerproduct/sobolev.h b/specific/model/innerproduct/sobolev.h index 1d759be..5f8ff1b 100644 --- a/specific/model/innerproduct/sobolev.h +++ b/specific/model/innerproduct/sobolev.h @@ -23,7 +23,7 @@ namespace sobolev double _bulkWeight; double _boundaryWeight; public: - 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.) {}; + InnerProduct(bool complex, const ConfigurationInterface* const config, const model::Discretization& m_discret, const gmshfem::common::GmshFem& gmshFem, std::string suffix=""); virtual void link(ModelUpdater* const mu, SensitivityUpdater* const mpu) override; virtual void unlink() override; virtual void modelIsObsolete() override; diff --git a/specific/optimization/descentsearch/gradientdescent.cpp b/specific/optimization/descentsearch/gradientdescent.cpp new file mode 100644 index 0000000..ba660e0 --- /dev/null +++ b/specific/optimization/descentsearch/gradientdescent.cpp @@ -0,0 +1,18 @@ +//GmshFWI Library +#include "gradientdescent.h" + +using namespace gmshfem; +using namespace gmshfem::common; + +namespace gradientdescent +{ + /* + * class DescentSearch + */ + void DescentSearch::operator()(FunctionalInterface* const functional) const + { + ModelField p = functional->gradient(); + p.oppose(); + functional->setModelPerturbation(p); + } +} diff --git a/specific/optimization/descentsearch/gradientdescent.h b/specific/optimization/descentsearch/gradientdescent.h new file mode 100644 index 0000000..f27fced --- /dev/null +++ b/specific/optimization/descentsearch/gradientdescent.h @@ -0,0 +1,40 @@ +#ifndef H_SPECIFIC_OPTIMIZATION_DESCENTSEARCH_GRADIENTSDESCENT +#define H_SPECIFIC_OPTIMIZATION_DESCENTSEARCH_GRADIENTSDESCENT + +//GmshFEM Library +#include "GmshFem.h" +#include "CSVio.h" + +//GmshFWI Library +#include "../../../common/optimization/descentsearch.h" + +namespace gradientdescent +{ + + /* + * class DescentSearchHistory + */ + class DescentSearchHistory final: public DescentSearchHistoryInterface + { + public: + virtual void write(std::string filename) const {}; + }; + + /* + * class DescentSearch + */ + class DescentSearch final: public DescentSearchInterface + { + private: + DescentSearchHistory* const _history; + public: + DescentSearch(const gmshfem::common::GmshFem& gmshFem) : _history( new DescentSearchHistory() ) {}; + ~DescentSearch() {delete _history;}; + + virtual void operator()(FunctionalInterface* const functional) const; + + virtual const DescentSearchHistoryInterface* const history() const {return _history;}; + }; +}; + +#endif // H_SPECIFIC_OPTIMIZATION_DESCENTSEARCH_GRADIENTSDESCENT diff --git a/specific/optimization/descentsearch/newDescentSearch.h b/specific/optimization/descentsearch/newDescentSearch.h new file mode 100644 index 0000000..bb789e3 --- /dev/null +++ b/specific/optimization/descentsearch/newDescentSearch.h @@ -0,0 +1,31 @@ +#ifndef H_SPECIFIC_OPTIMIZATION_NEWDESCENTSEARCH +#define H_SPECIFIC_OPTIMIZATION_NEWDESCENTSEARCH + +//Standard Library +#include <string> + +//GmshFEM Library +#include "GmshFem.h" +#include "Exception.h" + +//GmshFWI Library +#include "gradientdescent.h" +#include "newton_conjugategradient.h" + +DescentSearchInterface* newDescentSearch(const gmshfem::common::GmshFem& gmshFem) +{ + std::string descentsearch; + if(!gmshFem.userDefinedParameter(descentsearch, "descentsearch")) + { + throw gmshfem::common::Exception("Descent search type could not be found."); + } + if(descentsearch=="gradientdescent"){return new gradientdescent::DescentSearch(gmshFem);} + if(descentsearch=="newton_conjugategradient"){return new newton_conjugategradient::DescentSearch(gmshFem);} + else + { + throw gmshfem::common::Exception("Descent search " + descentsearch + " is not valid."); + return nullptr; + } +}; + +#endif //H_SPECIFIC_OPTIMIZATION_NEWDESCENTSEARCH diff --git a/specific/optimization/descentsearch/newton_conjugategradient.cpp b/specific/optimization/descentsearch/newton_conjugategradient.cpp new file mode 100644 index 0000000..bb27678 --- /dev/null +++ b/specific/optimization/descentsearch/newton_conjugategradient.cpp @@ -0,0 +1,139 @@ +//GmshFEM Library +#include "FieldRoutines.h" +//GmshFWI Library +#include "newton_conjugategradient.h" + +using namespace gmshfem; +using namespace gmshfem::common; + +static const std::complex< double > im = std::complex< double >(0., 1.); + +namespace newton_conjugategradient +{ + /* + * class DescentSearchHistoryLine2 + */ + void DescentSearchHistoryLine2::write(CSVio& file, bool write_header) const + { + if(write_header) + { + file << "Residual" << "Alpha" << "Beta" << csv::endl; + } + else + { + file << residual << alpha << beta << csv::endl; + } + } + + /* + * class DescentSearchHistoryLine + */ + void DescentSearchHistoryLine::write(CSVio& file, bool write_header) const + { + if(write_header) + { + file << "Success" << "Eta squared" << "Gradient norm squared" << csv::endl; + file << "History" << csv::endl; + } + else + { + file << success << eta2 << gradientNorm2 << csv::endl; + if(!history.empty()) + { + history.begin()->write(file,true); + for (auto it = history.begin(); it != history.end(); it++) + { + it->write(file); + } + } + } + } + + /* + * class DescentSearchHistory + */ + void DescentSearchHistory::write(std::string filename) const + { + CSVio file(filename, ';', common::OpeningMode::NewFile); + if(this->empty()) + { + file << "Descent 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 DescentSearch + */ + DescentSearch::DescentSearch(const GmshFem& gmshFem) : _history( new DescentSearchHistory() ) + { + if(!gmshFem.userDefinedParameter(_eta20, "eta20")) + { + throw Exception("Initial value of forcing sequence (eta20) could not be found."); + } + + if(!gmshFem.userDefinedParameter(_maxIteration, "descent_maxIteration")) + { + throw Exception("Maximum number of iterations for descent search (descent_maxIteration) could not be found."); + } + } + + void DescentSearch::operator()(FunctionalInterface* const functional) const + { + msg::print <<"Newton conjugate gradient" << msg::endl; + std::vector<DescentSearchHistoryLine2> historyline; + + double jj = std::real( functional->innerproduct()->product(functional->gradient(),functional->gradient()) ); + double eta2 = _eta20; + if (_history->size()>1) + { + double jj_1 = _history->back().gradientNorm2; + eta2 = std::max( _eta20, jj / jj_1 ); + } + msg::print << "\t" << "eta2 = " << eta2 << msg::endl; + msg::print << "\t" << "residual threshold = " << eta2 * jj << msg::endl; + + + ModelField x = functional->gradient(); + x.setValuesToZero(); + + ModelField r = functional->gradient(); + double rr = jj; + double rr_1 = rr; + + ModelField p = functional->gradient(); + p.setValuesToZero(); + + unsigned int n = 0; + while(rr > eta2 * jj && n < _maxIteration) + { + msg::print << "\t" << "--- Iteration #"+ std::to_string(n) +" ----" << msg::endl; + msg::print << "\t" << "residual = " << rr << msg::endl; + double beta = rr / rr_1; + multiplyAddAssignement(beta+ 0.*im, p, -1. + 0.*im,r); + msg::print << "\t \t" << "beta = " << beta << msg::endl; + + double alpha = rr / std::real( functional->innerproduct()->product(p,functional->hessian(p)) ); + msg::print << "\t \t" << "alpha = " << alpha << msg::endl; + multiplyAddAssignement(1.+0.*im, x, alpha + 0.*im,p); + multiplyAddAssignement(1.+0.*im, r, alpha + 0.*im,functional->hessian()); + + historyline.emplace_back(rr,alpha,beta); + + rr_1 = rr; + rr = std::real( functional->innerproduct()->product(r,r) ); + n++; + } + msg::print << "\t" << "residual = " << rr << msg::endl; + functional->setModelPerturbation(x); + _history->emplace_back(rr <= eta2 * jj, eta2, jj, historyline); + } +} diff --git a/specific/optimization/descentsearch/newton_conjugategradient.h b/specific/optimization/descentsearch/newton_conjugategradient.h new file mode 100644 index 0000000..7a9f6a0 --- /dev/null +++ b/specific/optimization/descentsearch/newton_conjugategradient.h @@ -0,0 +1,72 @@ +#ifndef H_SPECIFIC_OPTIMIZATION_DESCENTSEARCH_NEWTONCONJUGATEGRADIENT +#define H_SPECIFIC_OPTIMIZATION_DESCENTSEARCH_NEWTONCONJUGATEGRADIENT + +//GmshFEM Library +#include "GmshFem.h" +#include "CSVio.h" + +//GmshFWI Library +#include "../../../common/optimization/descentsearch.h" + +namespace newton_conjugategradient +{ + /* + * class DescentSearchHistoryLine2 + */ + class DescentSearchHistoryLine2 + { + public: + const double residual; + const double alpha; + const double beta; + + DescentSearchHistoryLine2(double residual_init, double alpha_init, double beta_init) : residual(residual_init), alpha(alpha_init), beta(beta_init) {}; + + void write(gmshfem::common::CSVio& file, bool write_header=false) const; + }; + + /* + * class DescentSearchHistoryLine + */ + class DescentSearchHistoryLine + { + public: + const bool success; + const double eta2; + const double gradientNorm2; + const std::vector<DescentSearchHistoryLine2> history; + + DescentSearchHistoryLine(bool success_init, double eta2_init, double gradientNorm2_init, const std::vector<DescentSearchHistoryLine2>& history_init) : success(success_init), eta2(eta2_init), gradientNorm2(gradientNorm2_init), history(history_init) {}; + + void write(gmshfem::common::CSVio& file, bool write_header=false) const; + }; + + /* + * class DescentSearchHistory + */ + class DescentSearchHistory final: public DescentSearchHistoryInterface, public std::vector<DescentSearchHistoryLine> + { + public: + virtual void write(std::string filename) const; + }; + + /* + * class DescentSearch + */ + class DescentSearch final: public DescentSearchInterface + { + private: + double _eta20; + unsigned int _maxIteration; + DescentSearchHistory* const _history; + public: + DescentSearch(const gmshfem::common::GmshFem& gmshFem); + ~DescentSearch() {delete _history;}; + + virtual void operator()(FunctionalInterface* const functional) const; + + virtual const DescentSearchHistoryInterface* const history() const {return _history;}; + }; +}; + +#endif // H_SPECIFIC_OPTIMIZATION_DESCENTSEARCH_NEWTONCONJUGATEGRADIENT diff --git a/specific/optimization/linesearch/backtracking.cpp b/specific/optimization/linesearch/backtracking.cpp new file mode 100644 index 0000000..61cd292 --- /dev/null +++ b/specific/optimization/linesearch/backtracking.cpp @@ -0,0 +1,100 @@ +//GmshFEM Library +#include "Exception.h" +#include "FieldRoutines.h" + +//GmshFWI Library +#include "backtracking.h" + +using namespace gmshfem; +using namespace gmshfem::common; + +static const std::complex< double > im = std::complex< double >(0., 1.); + +namespace backtracking +{ + /* + * class LineSearchHistoryLine + */ + void LineSearchHistoryLine::write(CSVio& file, bool write_header) const + { + if(write_header) + { + file << "Performance" << "a0" << "aN" << "N" << csv::endl; + } + else + { + file << performance << a0 << aN << N << csv::endl; + } + } + + /* + * class LineSearchHistory + */ + void LineSearchHistory::write(std::string filename) const + { + CSVio file(filename, ';', common::OpeningMode::NewFile); + if(this->empty()) + { + file << "Line search history is empty." << csv::endl; + } + else + { + this->begin()->write(file,true); + for (auto it = this->begin(); it != this->end(); it++) + { + it->write(file); + } + } + } + + /* + * class LineSearch + */ + LineSearch::LineSearch(const GmshFem& gmshFem) : _history( new LineSearchHistory() ) + { + if(!( + gmshFem.userDefinedParameter(_rho, "rho") && + gmshFem.userDefinedParameter(_c, "c") + )) + { + throw Exception("A backtracking linesearch parameter could not be found."); + } + } + + double LineSearch::operator()(FunctionalInterface* const functional) const + { + msg::print << "\t" << "--- Linesearch (backtracking) ----" << msg::endl; + ModelField m = functional->m(); + ModelField p = functional->dm(); + + double j = functional->performance(); + double dj = functional->directional1(p); + + double a0; + if (_history->size()>1) + { + double j_1 = _history->back().performance; + a0=std::min(1.,(j-j_1)/dj); + } + else + { + a0=1. /* *//*-0.125*jn/dj*/; + } + + double a=a0; + unsigned int i=0; + while (true) + { + double ja = functional->performance<ModelField>(multiplyAdd(1.+0.*im, m, a+0.*im,p)); + msg::print << "\t \t" << "j + c*a*dj = " << j + _c*a*dj << "\t" << " ja = " << ja << msg::endl; + if( ja < j + _c*a*dj ){break;} + else + { + a *= _rho; + } + i++; + } + _history->emplace_back(j,a0,a,i); + return a; + } +} diff --git a/specific/optimization/linesearch/backtracking.h b/specific/optimization/linesearch/backtracking.h new file mode 100644 index 0000000..d30a812 --- /dev/null +++ b/specific/optimization/linesearch/backtracking.h @@ -0,0 +1,57 @@ +#ifndef H_SPECIFIC_OPTIMIZATION_LINESEARCH_BACKTRACKING +#define H_SPECIFIC_OPTIMIZATION_LINESEARCH_BACKTRACKING + +//GmshFEM Library +#include "GmshFem.h" +#include "CSVio.h" + +//GmshFWI Library +#include "../../../common/optimization/linesearch.h" + +namespace backtracking +{ + /* + * class LineSearchHistoryLine + */ + class LineSearchHistoryLine + { + public: + const double performance; + const double a0; + const double aN; + const unsigned int N; + + LineSearchHistoryLine(double performance_init, double a0_init, double aN_init, unsigned int N_init) : performance(performance_init), a0(a0_init), aN(aN_init), N(N_init) {}; + + void write(gmshfem::common::CSVio& file, bool write_header=false) const; + }; + + /* + * class LineSearchHistory + */ + class LineSearchHistory final: public LineSearchHistoryInterface, public std::vector<LineSearchHistoryLine> + { + public: + virtual void write(std::string filename) const; + }; + + /* + * class LineSearch + */ + class LineSearch final: public LineSearchInterface + { + private: + double _rho; + double _c; + LineSearchHistory* const _history; + public: + LineSearch(const gmshfem::common::GmshFem& gmshFem); + ~LineSearch() {delete _history;}; + + virtual double operator()(FunctionalInterface* const functional) const; + + virtual const LineSearchHistoryInterface* const history() const {return _history;}; + }; +}; + +#endif // H_SPECIFIC_OPTIMIZATION_LINESEARCH_CLASSIC diff --git a/specific/optimization/linesearch/constant.cpp b/specific/optimization/linesearch/constant.cpp new file mode 100644 index 0000000..0f2f325 --- /dev/null +++ b/specific/optimization/linesearch/constant.cpp @@ -0,0 +1,73 @@ +//GmshFEM Library +#include "FieldRoutines.h" +#include "Exception.h" + +//GmshFWI Library +#include "constant.h" + +using namespace gmshfem; +using namespace gmshfem::common; + +static const std::complex< double > im = std::complex< double >(0., 1.); + +namespace constant +{ + /* + * class LineSearchHistoryLine + */ + void LineSearchHistoryLine::write(CSVio& file, bool write_header) const + { + if(write_header) + { + file << "Performance" << csv::endl; + } + else + { + file << performance << csv::endl; + } + } + + /* + * class LineSearchHistory + */ + void LineSearchHistory::write(std::string filename) const + { + CSVio file(filename, ';', common::OpeningMode::NewFile); + if(this->empty()) + { + file << "Line search history is empty." << csv::endl; + } + else + { + this->begin()->write(file,true); + for (auto it = this->begin(); it != this->end(); it++) + { + it->write(file); + } + } + } + + /* + * class LineSearch + */ + LineSearch::LineSearch(const GmshFem& gmshFem) : _history( new LineSearchHistory() ) + { + if(!( + gmshFem.userDefinedParameter(_a, "a") + )) + { + throw Exception("A constant linesearch parameter could not be found."); + } + } + + double LineSearch::operator()(FunctionalInterface* const functional) const + { + msg::print << "\t" << "--- Linesearch (constant) ----" << msg::endl; + ModelField m = functional->m(); + ModelField p = functional->dm(); + double j = functional->performance(); + functional->setModel( ((ModelField) multiplyAdd(1.+0.*im, m, _a+0.*im,p)) ); + _history->emplace_back(j); + return _a; + } +} diff --git a/specific/optimization/linesearch/constant.h b/specific/optimization/linesearch/constant.h new file mode 100644 index 0000000..24dbc62 --- /dev/null +++ b/specific/optimization/linesearch/constant.h @@ -0,0 +1,53 @@ +#ifndef H_SPECIFIC_OPTIMIZATION_LINESEARCH_CONSTANT +#define H_SPECIFIC_OPTIMIZATION_LINESEARCH_CONSTANT + +//GmshFEM Library +#include "GmshFem.h" +#include "CSVio.h" + +//GmshFWI Library +#include "../../../common/optimization/linesearch.h" + +namespace constant +{ + /* + * class LineSearchHistoryLine + */ + class LineSearchHistoryLine + { + public: + const double performance; + + LineSearchHistoryLine(double performance_init) : performance(performance_init) {}; + + void write(gmshfem::common::CSVio& file, bool write_header=false) const; + }; + + /* + * class LineSearchHistory + */ + class LineSearchHistory final: public LineSearchHistoryInterface, public std::vector<LineSearchHistoryLine> + { + public: + virtual void write(std::string filename) const; + }; + + /* + * class LineSearch + */ + class LineSearch final: public LineSearchInterface + { + private: + double _a; + LineSearchHistory* const _history; + public: + LineSearch(const gmshfem::common::GmshFem& gmshFem); + ~LineSearch() {delete _history;}; + + virtual double operator()(FunctionalInterface* const functional) const; + + virtual const LineSearchHistoryInterface* const history() const {return _history;}; + }; +}; + +#endif // H_SPECIFIC_OPTIMIZATION_LINESEARCH_CONSTANT diff --git a/specific/optimization/linesearch/newLineSearch.h b/specific/optimization/linesearch/newLineSearch.h new file mode 100644 index 0000000..50aea7c --- /dev/null +++ b/specific/optimization/linesearch/newLineSearch.h @@ -0,0 +1,33 @@ +#ifndef H_SPECIFIC_OPTIMIZATION_NEWLINESEARCH +#define H_SPECIFIC_OPTIMIZATION_NEWLINESEARCH + +//Standard Library +#include <string> + +//GmshFem Library +#include "GmshFem.h" +#include "Exception.h" + +//GmshFWI Library +#include "backtracking.h" +#include "secondorder.h" +#include "constant.h" + +LineSearchInterface* newLineSearch(const gmshfem::common::GmshFem& gmshFem) +{ + std::string linesearch; + if(!gmshFem.userDefinedParameter(linesearch, "linesearch")) + { + throw gmshfem::common::Exception("Line search type could not be found."); + } + if(linesearch=="backtracking"){return new backtracking::LineSearch(gmshFem);} + else if(linesearch=="secondorder"){return new secondorder::LineSearch(gmshFem);} + else if(linesearch=="constant"){return new constant::LineSearch(gmshFem);} + else + { + throw gmshfem::common::Exception("Line search " + linesearch + " is not valid."); + return nullptr; + } +}; + +#endif //H_SPECIFIC_OPTIMIZATION_NEWLINESEARCH diff --git a/specific/optimization/linesearch/secondorder.cpp b/specific/optimization/linesearch/secondorder.cpp new file mode 100644 index 0000000..949fa04 --- /dev/null +++ b/specific/optimization/linesearch/secondorder.cpp @@ -0,0 +1,72 @@ +//GmshFEM Library +#include "FieldRoutines.h" +#include "Exception.h" + +//GmshFWI Library +#include "secondorder.h" + +using namespace gmshfem; +using namespace gmshfem::common; + +static const std::complex< double > im = std::complex< double >(0., 1.); + +namespace secondorder +{ + /* + * class LineSearchHistoryLine + */ + void LineSearchHistoryLine::write(CSVio& file, bool write_header) const + { + if(write_header) + { + file << "Performance" << "a" << csv::endl; + } + else + { + file << performance << a << csv::endl; + } + } + + /* + * class LineSearchHistory + */ + void LineSearchHistory::write(std::string filename) const + { + CSVio file(filename, ';', common::OpeningMode::NewFile); + if(this->empty()) + { + file << "Line search history is empty." << csv::endl; + } + else + { + this->begin()->write(file,true); + for (auto it = this->begin(); it != this->end(); it++) + { + it->write(file); + } + } + } + + /* + * class LineSearch + */ + LineSearch::LineSearch(const GmshFem& gmshFem) : _history( new LineSearchHistory() ) + {} + + double LineSearch::operator()(FunctionalInterface* const functional) const + { + msg::print << "\t" << "--- Linesearch (second order) ----" << msg::endl; + ModelField m = functional->m(); + ModelField p = functional->dm(); + + double j = functional->performance(); + double dj = functional->directional1(p); + double d2j = functional->directional2(p,p); + double a = - dj / d2j; + + msg::print << "\t \t" << "a = " << a << msg::endl; + functional->setModel( ((ModelField) multiplyAdd(1.+0.*im, m, a+0.*im,p)) ); + _history->emplace_back(j,a); + return a; + } +} diff --git a/specific/optimization/linesearch/secondorder.h b/specific/optimization/linesearch/secondorder.h new file mode 100644 index 0000000..2696106 --- /dev/null +++ b/specific/optimization/linesearch/secondorder.h @@ -0,0 +1,53 @@ +#ifndef H_SPECIFIC_OPTIMIZATION_LINESEARCH_SECONDORDER +#define H_SPECIFIC_OPTIMIZATION_LINESEARCH_SECONDORDER + +//GmshFEM Library +#include "GmshFem.h" +#include "CSVio.h" + +//GmshFWI Library +#include "../../../common/optimization/linesearch.h" + +namespace secondorder +{ + /* + * class LineSearchHistoryLine + */ + class LineSearchHistoryLine + { + public: + const double performance; + const double a; + + LineSearchHistoryLine(double performance_init, double a_init) : performance(performance_init), a(a_init) {}; + + void write(gmshfem::common::CSVio& file, bool write_header=false) const; + }; + + /* + * class LineSearchHistory + */ + class LineSearchHistory final: public LineSearchHistoryInterface, public std::vector<LineSearchHistoryLine> + { + public: + virtual void write(std::string filename) const; + }; + + /* + * class LineSearch + */ + class LineSearch final: public LineSearchInterface + { + private: + LineSearchHistory* const _history; + public: + LineSearch(const gmshfem::common::GmshFem& gmshFem); + ~LineSearch() {delete _history;}; + + virtual double operator()(FunctionalInterface* const functional) const; + + virtual const LineSearchHistoryInterface* const history() const {return _history;}; + }; +}; + +#endif // H_SPECIFIC_OPTIMIZATION_LINESEARCH_SECONDORDER diff --git a/specific/optimization/minimumsearch/classic.cpp b/specific/optimization/minimumsearch/classic.cpp index ce777a0..f1ce4ed 100644 --- a/specific/optimization/minimumsearch/classic.cpp +++ b/specific/optimization/minimumsearch/classic.cpp @@ -76,13 +76,13 @@ namespace classic */ MinimumSearch::MinimumSearch(const GmshFem& gmshFem) : _history( new MinimumSearchHistory() ) { - if(!gmshFem.userDefinedParameter(_maxIteration, "maxIteration")) + if(!gmshFem.userDefinedParameter(_maxIteration, "minimum_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 + void MinimumSearch::operator()(ModelField* const m, FunctionalInterface* const functional, const DescentSearchInterface* const descentsearch, const LineSearchInterface* const linesearch) const { std::vector<MinimumSearchHistoryLine2> historyline; @@ -95,7 +95,7 @@ namespace classic 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)); + //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); @@ -103,9 +103,10 @@ namespace classic if(n >= _maxIteration){break;} // ou critère sur djn ou critère sur jn else { - //descent_direction(); - //linesearch(); - //m = m+a*p + (*descentsearch)(functional); + (*linesearch)(functional); + *m = functional->m(); + //m->write("m"+std::to_string(n)); } n++; } diff --git a/specific/optimization/minimumsearch/classic.h b/specific/optimization/minimumsearch/classic.h index eb04634..40bb363 100644 --- a/specific/optimization/minimumsearch/classic.h +++ b/specific/optimization/minimumsearch/classic.h @@ -60,7 +60,7 @@ namespace classic 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 void operator()(ModelField* const m, FunctionalInterface* const functional, const DescentSearchInterface* const descentsearch, const LineSearchInterface* const linesearch) const; virtual const MinimumSearchHistoryInterface* const history() const {return _history;}; }; diff --git a/specific/wave/discretization.cpp b/specific/wave/discretization.cpp index 8f5f23b..9b935f2 100644 --- a/specific/wave/discretization.cpp +++ b/specific/wave/discretization.cpp @@ -14,14 +14,23 @@ namespace wave /* * class Discretization */ - Discretization<Physic::acoustic>::Discretization(const GmshFem& gmshFem) + Discretization<Physic::acoustic>::Discretization(const GmshFem& gmshFem, std::string suffix) { msg::print << "Read wave discretization parameters" << msg::endl; std::string str_buffer; - if( - !(gmshFem.userDefinedParameter(_functionSpaceDegree, "wave_FunctionSpaceDegree") && - gmshFem.userDefinedParameter(str_buffer, "wave_FunctionSpaceType")) - ) + if(!( + ( + gmshFem.userDefinedParameter(_functionSpaceDegree, "wave_FunctionSpaceDegree") + || + gmshFem.userDefinedParameter(_functionSpaceDegree, "wave_FunctionSpaceDegree"+suffix) + ) + && + ( + gmshFem.userDefinedParameter(str_buffer, "wave_FunctionSpaceType") + || + gmshFem.userDefinedParameter(str_buffer, "wave_FunctionSpaceType"+suffix) + ) + )) { throw common::Exception("A discretization parameter could not be found."); } diff --git a/specific/wave/discretization.h b/specific/wave/discretization.h index 7d4d4be..8c7e059 100644 --- a/specific/wave/discretization.h +++ b/specific/wave/discretization.h @@ -23,7 +23,7 @@ namespace wave 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;}; diff --git a/specific/wave/equation/helmholtz.cpp b/specific/wave/equation/helmholtz.cpp index 5bed9bc..4887459 100644 --- a/specific/wave/equation/helmholtz.cpp +++ b/specific/wave/equation/helmholtz.cpp @@ -85,9 +85,11 @@ namespace helmholtz { _formulation.galerkin(-d.state(Type::PA).value(s,r), tf(_v), _config->receiver(s,r), integrationType(0)); } - _formulation.galerkin( lambda[s] * _pulsation * _pulsation * dm, tf(_v), _config->model_unknown(Support::BLK), integrationType(2*degree)); - - _formulation.galerkin(-im * lambda[s] * _pulsation / 2. / sqrt(m) * dm, tf(_v), _config->model_unknown(Support::BND), integrationType(2*degree)); //-1 from sign conventions + if(!_gaussNewton) + { + _formulation.galerkin( lambda[s] * _pulsation * _pulsation * dm, tf(_v), _config->model_unknown(Support::BLK), integrationType(2*degree)); + _formulation.galerkin(-im * lambda[s] * _pulsation / 2. / sqrt(m) * dm, tf(_v), _config->model_unknown(Support::BND), integrationType(2*degree)); //-1 from sign conventions + } break; } } @@ -113,18 +115,38 @@ namespace helmholtz switch (support) { case Support::BLK: - return - _pulsation * _pulsation * - ( - correlate(u,dlambda) - +correlate(du,lambda) - ); + if(!_gaussNewton) + { + return - _pulsation * _pulsation * + ( + correlate(u,dlambda) + +correlate(du,lambda) + ); + } + else + { + return - _pulsation * _pulsation * + ( + correlate(u,dlambda) + ); + } case Support::BND: - return im * _pulsation / 2. / sqrt(m) * - ( - correlate(u,dlambda) - +correlate(du,lambda) - -dm / 2. / m * correlate(u,lambda) - ); + if(!_gaussNewton) + { + return im * _pulsation / 2. / sqrt(m) * + ( + correlate(u,dlambda) + +correlate(du,lambda) + -dm / 2. / m * correlate(u,lambda) + ); + } + else + { + return im * _pulsation / 2. / sqrt(m) * + ( + correlate(u,dlambda) + ); + } } break; case Order::DIAG: @@ -146,7 +168,9 @@ namespace helmholtz updateGreenReceiver(ms); double pulsation4 = _pulsation * _pulsation * _pulsation * _pulsation; double wl = wavelength(_config->m_reference()); - double integrationSupport = M_PI * (wl/4.) * (wl/4.); + double integrationSupport = M_PI * (wl/2.) * (wl/2.); + //wl/4. signifie que les interférences destructive apparaissent en xr->+infty et xr -> -infty (ou à l'opposé si les récepteur sont sur un cercle) + //wl/2. signifie que la phase du préconditionneur fait un tour entier de xr->+infty à xr -> -infty (ou un tour si les récepteur sont sur un cercle) return pulsation4 * integrationSupport * green0_preconditioner<Physic::acoustic>(_g,_config); } }; diff --git a/specific/wave/equation/helmholtz.h b/specific/wave/equation/helmholtz.h index 0ec2e0d..013f915 100644 --- a/specific/wave/equation/helmholtz.h +++ b/specific/wave/equation/helmholtz.h @@ -15,8 +15,14 @@ namespace helmholtz { private: double _pulsation; + bool _gaussNewton; public: - Equation(double pulsation,const ConfigurationInterface* const config, const wave::Discretization<Physic::acoustic>& w_discret,const gmshfem::common::GmshFem& gmshFem) : DifferentialEquationInterface<Physic::acoustic>(config,w_discret,gmshFem), _pulsation(pulsation) {}; + Equation(double pulsation,const ConfigurationInterface* const config, const wave::Discretization<Physic::acoustic>& w_discret,const gmshfem::common::GmshFem& gmshFem) : DifferentialEquationInterface<Physic::acoustic>(config,w_discret,gmshFem), _pulsation(pulsation) + { + unsigned int gn = 0; + gmshFem.userDefinedParameter(gn, "gaussNewton"); + _gaussNewton = ((bool) gn); + }; //Model related virtual ModelFunction update_sensitivity(Order order, Support support, const ModelStateEvaluator& ms, const WaveStateEvaluator<Physic::acoustic>& ws); diff --git a/synthetics.cpp b/synthetics.cpp index e3c121b..db00d8d 100644 --- a/synthetics.cpp +++ b/synthetics.cpp @@ -17,27 +17,36 @@ int synthetics(const GmshFem& gmshFem) { std::string name = "noname"; gmshFem.userDefinedParameter(name, "name"); - ConfigurationInterface* configuration = newConfiguration(gmshFem); + ConfigurationInterface* configuration = newConfiguration(name, gmshFem); save(configuration->m0(), configuration->model_known(Support::BLK), name+"_m_syn_blk", "pos", ""); save(configuration->m0(), configuration->model_known(Support::BND), name+"_m_syn_bnd", "pos", ""); - wave::Discretization<T_Physic> w_discret(gmshFem); - double frequency = 0.; - if(!gmshFem.userDefinedParameter(frequency, "frequency")) - { - throw Exception("Frequency could not be found."); - } - msg::print << "Generate synthetic data at frequency = " << frequency << msg::endl; - EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem); - WaveUpdater<T_Physic> wu(configuration,nullptr,equation); - DataUpdater<T_Physic> du(configuration,&wu,nullptr); ModelState ms(configuration); - std::array<bool,4> needUpToDate = {true,false,false,false}; - ((du.get(needUpToDate,ms)).state(Type::F)).write(name+"_data"); - wu.get(needUpToDate,ms).write(Type::F,name+"_u_syn"); + unsigned int n_freq = 1; + gmshFem.userDefinedParameter(n_freq, "n_freq"); + for (unsigned int f = 0; f < n_freq; f++) + { + std::string suffix = std::to_string(f); + + double frequency; + if(!gmshFem.userDefinedParameter(frequency, "frequency"+suffix)) + { + throw Exception("Frequency "+suffix+" could not be found."); + } + msg::print << "Generate synthetic data at frequency #"+suffix+" = " << frequency << msg::endl; - delete equation; + wave::Discretization<T_Physic> w_discret(gmshFem,suffix); + EquationInterface<T_Physic>* const equation = newEquation(2.*M_PI*frequency,configuration,w_discret,gmshFem); + WaveUpdater<T_Physic> wu(configuration,nullptr,equation); + DataUpdater<T_Physic> du(configuration,&wu,nullptr); + + std::array<bool,4> needUpToDate = {true,false,false,false}; + ((du.get(needUpToDate,ms)).state(Type::F)).write(name+"_data"+suffix); + wu.get(needUpToDate,ms).write(Type::F,name+"_u_syn"+suffix); + + delete equation; + } delete configuration; return 0; -- GitLab