Flooding regularization

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

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

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));
+plt.xlabel('receiver index [-]')
+plt.ylabel(r'time [s]')
\ No newline at end of file

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", "");

specific/model/regularization/flooding.cpp

+//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;
+            }
+        }
+    }
+};

specific/model/regularization/flooding.h

+#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

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.");