Select Git revision
simplewave.cpp
-
Boris Martin authoredBoris Martin authored
simplewave.cpp 12.01 KiB
//Standard library
//GmshFem library
#include "Message.h"
#include "Exception.h"
#include "Function.h"
#include "Post.h"
//GmshFWI Library
#include "simplewave.h"
#include "../../model/parametrization/simplewave.h"
#include "../../configuration/green0_preconditioner.h"
#include "../correlation.h"
using namespace gmshfem;
using namespace gmshfem::common;
using namespace gmshfem::equation;
using namespace gmshfem::function;
using namespace gmshfem::post;
static const std::complex< double > im = std::complex< double >(0., 1.);
/*
* Simple wave (acoustics)
*/
#define u ws.state(Type::FS)
#define du ws.state(Type::PFS)
#define lambda ws.state(Type::AS)
#define dlambda ws.state(Type::PAS)
#define m0 _config->m0_natural(_pulsation)[0]
#define m ms.state(Type::FS)[0]
#define dm ms.state(Type::PFS)[0]
#define _w_degree _v.degree()
namespace simplewave
{
Equation::Equation(unsigned int f_idx,double pulsation,const ConfigurationInterface* const config, const wave::Discretization<Physic::acoustic>& w_discret,const gmshfem::common::GmshFem& gmshFem, std::string suffix) : DifferentialEquationInterface<Physic::acoustic>(config,w_discret,gmshFem,suffix), _f_idx(f_idx), _pulsation(pulsation)
{
unsigned int gn = 0;
if(!gmshFem.userDefinedParameter(gn, "gaussNewton"))
{
msg::warning << "Gauss-Newton approximation (equation) flag could not be found. Approximation is not used (default)." << msg::endl;
}
_gaussNewton = ((bool) gn);
unsigned int ri = 0;
if(!gmshFem.userDefinedParameter(ri, "equation_boundary_realImpedance"))
{
msg::warning << "Real boundary impedance flag could not be found. Complex boundary impedances are used (default)." << msg::endl;
}
_realImpedance = ((bool) ri);
unsigned int gp = 1;
if(!gmshFem.userDefinedParameter(gp, "equation_greenPreconditioner"))
{
msg::warning << "Green preconditionner (equation) flag could not be found. Green preconditioner is used (default)." << msg::endl;
}
_greenPreconditioner = ((bool) gp);
if(!_greenPreconditioner)
{
if(!gmshFem.userDefinedParameter(_const_preconditioner, "equation_const_preconditioner"))
{
throw Exception("Impossible to find constant preconditioner (equation)");
}
}
if(_realImpedance && _boundary)
{
throw Exception("Impossible to compute boundary terms with real impedance: sensitivities are wrong.");
}
if (gmshFem.userDefinedParameter(_topBC, "top_bc")) {
if (_topBC != "Absorbing" && _topBC != "Dirichlet" && _topBC != "Neumann") {
throw Exception("Invalid value of parameter top_bc: " + _topBC);
}
}
else {
msg::warning << "SimpleWave without specified top_bc defaults to Absorbing." << msg::endl;
}
}
bool Equation::modelIsObsolete()
{
_systemIsUpToDate=false;
_systemIsFactorized=false;
_geAreUpToDate=false;
_grAreUpToDate=false;
return _greenPreconditioner;
}
bool Equation::compatible(const ParametrizationInterface* const parametrization) const
{
return
(
nullptr != dynamic_cast<const natural::Parametrization* const>(parametrization)
||
nullptr != dynamic_cast<const slowness2_real::Parametrization* const>(parametrization)
||
nullptr != dynamic_cast<const ln_slowness2_real::Parametrization* const>(parametrization)
||
nullptr != dynamic_cast<const permcond_bireal::Parametrization* const>(parametrization)
||
nullptr != dynamic_cast<const ln_permcond_bireal::Parametrization* const>(parametrization)
);
}
double Equation::wavelength(std::vector<std::complex<double>> mc) const
{
return 2. * M_PI / (_pulsation * std::real(std::sqrt(mc[0])));
}
void Equation::setLHS(const ModelStateEvaluator& ms)
{
auto top = _config->named_domains().find("top_bnd");
if (_topBC == "Dirichlet" || _topBC == "Neumann")
{
if (top == _config->named_domains().end())
throw Exception("Asked for specific BCs on top but this entity isn't defined.");
}
if (_topBC == "Dirichlet")
{
_v.addConstraint(top->second, 0.0);
}
_formulation.integral(-1.*grad(dof(_v)), grad(tf(_v)), _config->wave_omega(Support::BLK),integrationType(2*_w_degree));
_formulation.integral(_pulsation * _pulsation * m0 * dof(_v), tf(_v), _config->model_known(Support::BLK), integrationType(_integrationDegreeBlk));
auto absorbingPartKnown = _config->model_known(Support::BND);
auto absorbingPartUnknown = _config->model_unknown(Support::BND);
if (_topBC == "Neumann") {
absorbingPartKnown &= ~(top->second);
absorbingPartUnknown &= ~(top->second);
}
if(_realImpedance)
{
_formulation.integral(-1. * im * _pulsation * sqrt(real(m0)) * dof(_v), tf(_v), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions }
else
{
_formulation.integral(-1. * im * _pulsation * sqrt(m0) * dof(_v), tf(_v), absorbingPartKnown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
if(!_config->model_unknown(Support::BLK).isEmpty())
{
_formulation.integral(_pulsation * _pulsation * m * dof(_v), tf(_v), _config->model_unknown(Support::BLK), integrationType(_integrationDegreeBlk));
}
if(!_config->model_unknown(Support::BND).isEmpty())
{
if(_realImpedance)
{
_formulation.integral(-1. * im * _pulsation * sqrt(real(m)) * dof(_v), tf(_v), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
else
{
_formulation.integral(-1. * im * _pulsation * sqrt(m) * dof(_v), tf(_v), absorbingPartUnknown, integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
}
}
void Equation::setRHS(Type type, unsigned int s, const DataStateEvaluator<Physic::acoustic>& d, const ModelStateEvaluator& ms, const WaveStateEvaluator<Physic::acoustic>& ws)
{
switch (type)
{
case Type::FS:
for (unsigned int e = 0; e < _config->ne(s); e++)
{
_formulation.integral(-1.,tf(_v), _config->emitter(s,e), integrationType(0));
}
break;
case Type::AS:
for (unsigned int r = 0; r < _config->nr(s); r++)
{
_formulation.integral(-d.state(Type::AS).value(_f_idx,s,r),tf(_v), _config->receiver(s,r), integrationType(0));
}
break;
case Type::PFS:
_formulation.integral( u[s] * _pulsation * _pulsation * dm, tf(_v), _config->model_unknown(Support::BLK), integrationType(_integrationDegreeBlk));
if(_boundary)
{
_formulation.integral(-im * u[s] * _pulsation / 2. / sqrt(m)* dm, tf(_v), _config->model_unknown(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
break;
case Type::PAS:
for (unsigned int r = 0; r < _config->nr(s); r++)
{
_formulation.integral(-d.state(Type::PAS).value(_f_idx,s,r), tf(_v), _config->receiver(s,r), integrationType(0));
}
if(!_gaussNewton)
{
_formulation.integral( lambda[s] * _pulsation * _pulsation * dm, tf(_v), _config->model_unknown(Support::BLK), integrationType(_integrationDegreeBlk));
if(_boundary)
{
_formulation.integral(-im * lambda[s] * _pulsation / 2. / sqrt(m) * dm, tf(_v), _config->model_unknown(Support::BND), integrationType(_integrationDegreeBnd)); //-1 from sign conventions
}
}
break;
case Type::DS: default:
throw Exception("Impossible to set RHS other than F,PF,A,PA");
break;
}
}
void Equation::setGreenRHS(unsigned int p)
{
_formulation.integral(-1.,tf(_v), _config->point(p), integrationType(0));
}
Sensitivity Equation::update_sensitivity(Order order, Support support, const DataStateEvaluator<Physic::acoustic>& ds, const ModelStateEvaluator& ms, const WaveStateEvaluator<Physic::acoustic>& ws)
{ ModelMonoFunction fct;
switch (order)
{
case Order::FST:
switch (support)
{
case Support::BLK:
fct = - _pulsation * _pulsation * correlate(u,lambda);
break;
case Support::BND:
fct = im * _pulsation / 2. / sqrt(m) * correlate(u,lambda);
break;
}
break;
case Order::SCD:
switch (support)
{
case Support::BLK:
if(!_gaussNewton)
{
fct = - _pulsation * _pulsation *
(
correlate(u,dlambda)
+correlate(du,lambda)
);
}
else
{
fct = - _pulsation * _pulsation *
(
correlate(u,dlambda)
);
}
break;
case Support::BND:
if(!_gaussNewton)
{
fct = im * _pulsation / 2. / sqrt(m) *
(
correlate(u,dlambda)
+correlate(du,lambda)
-dm / 2. / m * correlate(u,lambda)
);
}
else
{
fct = im * _pulsation / 2. / sqrt(m) *
(
correlate(u,dlambda)
);
}
break;
}
break;
case Order::DIAG:
switch (support)
{
case Support::BLK:
return preconditioner(ds,ms);
case Support::BND:
break;
}
break;
}
return Sensitivity(1,ModelMonoFunction(fct));
}
Sensitivity Equation::preconditioner(const DataStateEvaluator<Physic::acoustic>& ds, const ModelStateEvaluator& ms)
{
if(_greenPreconditioner) {
if(!_geAreUpToDate){updateGreenEmitter(ms);};
if(!_grAreUpToDate){updateGreenReceiver(ms);};
double pulsation4 = _pulsation * _pulsation * _pulsation * _pulsation;
double wl = wavelength(_config->mc_natural(_pulsation));
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 Sensitivity(1,green0_preconditioner<Physic::acoustic>(ds.state(Type::DS),_g,_config) * pulsation4 * integrationSupport);
}
else
{
return Sensitivity(1,_const_preconditioner);
}
}
};