diff --git a/input/VO2009Algorithm/acoustic/natural/common.txt b/input/VO2009Algorithm/acoustic/natural/common.txt index 4787d815894eec5d8770351f31deeabcced06264..197ffafec4a4f77d2805e09b4fdf9ec97cbc996b 100644 --- a/input/VO2009Algorithm/acoustic/natural/common.txt +++ b/input/VO2009Algorithm/acoustic/natural/common.txt @@ -6,7 +6,6 @@ equation_boundary = 0 # #Parametrization parametrization = slowness2_real -model_reference_frequency = 1. # #Configuration configuration = volume_acquisition diff --git a/input/VO2009Algorithm/acoustic/natural/synthetics.txt b/input/VO2009Algorithm/acoustic/natural/synthetics.txt index 8ac504a4db942bd4e389057890f79e94d3228221..bfd712418b4d4b7c1fe067dceaf921c947478b0d 100644 --- a/input/VO2009Algorithm/acoustic/natural/synthetics.txt +++ b/input/VO2009Algorithm/acoustic/natural/synthetics.txt @@ -14,5 +14,5 @@ Im(mi0c0) = 0. #Discretization h = 0.05 #Output -write_wave_fields = 0 +write_wave_fields = 1 write_data_fields = 1 diff --git a/plot_timewave.py b/plot_timewave.py index 52e944e4cf7fe2eba17b75f71bb6d811d4befe5b..1f20935ceaf76121af3dcf705179b062489a01bd 100644 --- a/plot_timewave.py +++ b/plot_timewave.py @@ -8,6 +8,8 @@ Created on Thu Aug 25 14:43:28 2022 import numpy as np import matplotlib.pyplot as plt from matplotlib import rc +plt.rcParams['ytick.right'] = plt.rcParams['ytick.labelright'] = False +plt.rcParams['ytick.left'] = plt.rcParams['ytick.labelleft'] = True plt.rcParams['xtick.bottom'] = plt.rcParams['xtick.labelbottom'] = False plt.rcParams['xtick.top'] = plt.rcParams['xtick.labeltop'] = True rc('text.latex', preamble=r'\usepackage{amsmath} \usepackage{amssymb} \renewcommand{\familydefault}{\sfdefault} \usepackage{sansmathfonts}') @@ -17,14 +19,13 @@ lw = 1.5 H = 2.0; L = 3.0; -""" # Marmousi path = "extended_marmousi_data/" -prename = "timewave_synthetics_data"; s = 61 -#prename = "timewave_synthetics_data_snr2"; s = 61 +#prename = "timewave_synthetics_data"; s = 61 +prename = "timewave_synthetics_data_snr2"; s = 61 #prename = "paper2_extended_timewave_initial_data"; s = 0 -prename = "paper2_extended_timewave_inverse_data"; s = 0 +#prename = "paper2_extended_timewave_inverse_data"; s = 0 T=8.0 Δt=0.004 @@ -36,7 +37,6 @@ Te=4.0 time_samp = 4 ΔNr = 60 ΔT = 1. -""" # Concrete 3 """ @@ -59,6 +59,7 @@ time_samp = 4 ΔT = 25 #[ms] """ +""" # Cross path = "cross_data/" @@ -78,7 +79,7 @@ Te=50 time_samp = 20 ΔNr = 20 ΔT = 10 - +""" """ # S-Cube @@ -110,9 +111,16 @@ for n in range(0,N+1,time_samp): data[np.int(n/time_samp),:] = tmp[s,0:2*Nr:2] figure=plt.figure(figsize=(L,2*H),tight_layout=False) -plt.subplots_adjust(top=0.92,right=0.9,bottom=0.03,left=0.1) +plt.subplots_adjust(top=0.93,right=0.9,bottom=0.06,left=0.1) +ax=plt.subplot() +ax.yaxis.set_label_position("right") +ax.set_ylabel('Y', rotation=-90, labelpad=15) +ax.xaxis.set_label_position("bottom") +ax.set_xlabel('Y', rotation=0, labelpad=5) plt.gca().invert_yaxis() time = np.arange(0,N+time_samp+1,time_samp)*Δt #[s] plt.pcolormesh(np.arange(0,Nr+1,1)-0.5,time-Δt*time_samp/2,data[:,:],vmin=-amplitude,vmax=amplitude,cmap='binary') plt.xticks(np.arange(0,Nr,ΔNr)); -plt.yticks(np.arange(0,Te*1.1,ΔT)); \ No newline at end of file +plt.yticks(np.arange(0,Te*1.1,ΔT)); +plt.xlabel('receiver index [-]') +plt.ylabel(r'time [s]') \ No newline at end of file diff --git a/specific/model/innerproduct/sobolev.cpp b/specific/model/innerproduct/sobolev.cpp index f157eb4d382de5d0096dab9dc78f2e95b473cd59..13f09da4e8a3850af54bf1c9ec04625be954cd4d 100644 --- a/specific/model/innerproduct/sobolev.cpp +++ b/specific/model/innerproduct/sobolev.cpp @@ -346,6 +346,21 @@ namespace sobolev_df ) / gradm2; } + /* + bool _flooding = true; + if(_flooding) + { + ScalarFunction< std::complex< double > > partition = heaviside< std::complex< double > >( yComp(grad(m[0]))+1e-7 ); + _D = partition * _D + + + (1-partition)*tensor< std::complex< double > > + ( + 0.,0.,0., + 0.,1.,0., + 0.,0.,0. + ); + } + */ /* ScalarFunction< std::complex< double > > gradm2 = ScalarFunction< std::complex< double > >(_m2x[0])+ScalarFunction< std::complex< double > >(_m2y[0]); gmshfem::post::save(sqrt(gradm2), _config->model_unknown(Support::BLK), "absgradm0", "pos", ""); diff --git a/specific/model/regularization/flooding.cpp b/specific/model/regularization/flooding.cpp new file mode 100644 index 0000000000000000000000000000000000000000..a47afd8616db0b1c08ad93333ad1bd251c2ae063 --- /dev/null +++ b/specific/model/regularization/flooding.cpp @@ -0,0 +1,82 @@ +//GmshFEM Library +#include "Function.h" + +//GmshFWI Library +#include "flooding.h" + +using namespace gmshfem; +using namespace gmshfem::common; +using namespace gmshfem::equation; +using namespace gmshfem::function; +using namespace gmshfem::post; + +/* +* Flooding +*/ +namespace flooding +{ + Regularization::Regularization(const ConfigurationInterface* const config, const ModelField& m, const gmshfem::common::GmshFem& gmshFem, std::string suffix) : DifferentialRegularizationInterface(config,m,gmshFem,suffix), _lambda(_m_size,0.) + { + for (unsigned int c = 0; c < _m_size; c++) + { + std::string suffix1 = "c"+std::to_string(c); + if(!( + gmshFem.userDefinedParameter(_lambda[c], "regularization_lambda") || + gmshFem.userDefinedParameter(_lambda[c], "regularization_lambda"+suffix) + || + gmshFem.userDefinedParameter(_lambda[c], "regularization_lambda"+suffix1) || + gmshFem.userDefinedParameter(_lambda[c], "regularization_lambda"+suffix1+suffix) + )) + { + throw Exception("Regularization scale could not be found."); + } + } + if(!( + gmshFem.userDefinedParameter(_p, "regularization_p") || + gmshFem.userDefinedParameter(_p, "regularization_p"+suffix) + )) + { + throw Exception("Regularization order could not be found."); + } + }; + + double Regularization::performance(const ModelField& m) + { + std::complex<double> out = 0.; + for (unsigned int c = 0; c < _m_size; c++) + { + ScalarFunction<std::complex<double>> mz = yComp(grad(m[c])); + ScalarFunction<std::complex<double>> f = (abs(mz) - mz) / 2.; + out += _lambda[c] * integrate( pow(f,2*_p) / 2. / _p, _config->model_unknown(Support::BLK), integrationType(_integrationDegreeBlk) ); + } + return std::real(out); + } + + void Regularization::setRHS(Order order, Support support, const ModelState& ms) + { + for (unsigned int c = 0; c < _m_size; c++) + { + ScalarFunction<std::complex<double>> mz = yComp(grad(ms.field(Type::FS)[c])); + ScalarFunction<std::complex<double>> f = (abs(mz) - mz) / 2.; + + ScalarFunction<std::complex<double>> partition = heaviside< std::complex< double > >( -mz ); + VectorFunction< std::complex< double > > integrand = vector< std::complex< double > >(0.,-partition * pow(f,2*_p-1),0.); + + switch (support) + { + case Support::BND:default:break; + case Support::BLK: + switch (order) + { + case Order::FST: + _formulation[c].integral(-1. * _lambda[c] * integrand, grad(tf(_k[c])),_config->model_unknown(Support::BLK), integrationType(_integrationDegreeBlk)); + break; + case Order::SCD: + break; + default:break; + } + break; + } + } + } +}; diff --git a/specific/model/regularization/flooding.h b/specific/model/regularization/flooding.h new file mode 100644 index 0000000000000000000000000000000000000000..fb9fd5d8ce6d7a86707eb2914a1043d1517a2fc5 --- /dev/null +++ b/specific/model/regularization/flooding.h @@ -0,0 +1,29 @@ +#ifndef H_SPECIFIC_WAVE_REGULARIZATION_FLOODING +#define H_SPECIFIC_WAVE_REGULARIZATION_FLOODING + +//GmshFEM Library +#include "GmshFem.h" + +//GmshFWI Library +#include "../../../common/model/regularization/regularization.h" + +/* +* Flooding +*/ +namespace flooding +{ + class Regularization final: public DifferentialRegularizationInterface + { + private: + double _p; + std::vector<double> _lambda; + public: + Regularization(const ConfigurationInterface* const config, const ModelField& m, const gmshfem::common::GmshFem& gmshFem, std::string suffix=""); + + virtual double performance(const ModelField& m); + + virtual void setRHS(Order order, Support support, const ModelState& ms); + }; +}; + +#endif //H_SPECIFIC_WAVE_REGULARIZATION_FLOODING diff --git a/specific/model/regularization/newRegularization.cpp b/specific/model/regularization/newRegularization.cpp index 28574b34341cb9abbaf4f21240fd42988252b692..1d259315584d86e05c55869260ba216bd613dc79 100644 --- a/specific/model/regularization/newRegularization.cpp +++ b/specific/model/regularization/newRegularization.cpp @@ -7,6 +7,7 @@ //GmshFWI Library #include "tikhonov.h" #include "totalvariation.h" +#include "flooding.h" RegularizationInterface* newRegularization(const ConfigurationInterface* const config, const ModelField& m, const gmshfem::common::GmshFem& gmshFem, std::string suffix) { @@ -18,6 +19,7 @@ RegularizationInterface* newRegularization(const ConfigurationInterface* const c } if(regularization=="tikhonov"){return new tikhonov::Regularization(config,m,gmshFem, suffix);} else if(regularization=="totalvariation"){return new totalvariation::Regularization(config,m,gmshFem, suffix);} + else if(regularization=="flooding"){return new flooding::Regularization(config,m,gmshFem, suffix);} else { throw gmshfem::common::Exception("Regularization " + regularization + " is not valid.");