diff --git a/FWI/src/FWI_DGTS.cpp b/FWI/src/FWI_DGTS.cpp
index 517769c0080c23417b08936ad640db5af66ef6be..911a45da3fe2a0428e97ada9cc1f39a6a39278c6 100644
--- a/FWI/src/FWI_DGTS.cpp
+++ b/FWI/src/FWI_DGTS.cpp
@@ -143,7 +143,7 @@ void createMeasuredSeismogram(const std::shared_ptr<ForwardModel>& FWD,WPConf &C
Config::GetDataPath() + std::string("FWI/Obs") + std::to_string(i));
}
}
- if (signalToNoise > -100000 + Eps) {
+ if (signalToNoise > -100000 + ForwardModel::Eps) {
int seisnr = 0;
for (auto &seis : obs) {
seis.getSeismogramWrite().addRelativeNoise(signalToNoise,
@@ -246,11 +246,11 @@ void fillInitialModel(ForwardModel &FWD, FWIDataContainerDomain &initial) {
for (cell c = M.cells(); c != M.cells_end(); ++c) {
double tmp0 = Tmp(c(), 0);
double tmp1 = Tmp(c(), 1);
- M(c(), 0) = tmp0 > Eps ? log(tmp0) : 0;
- M(c(), 1) = tmp1 > Eps ? log(tmp1) : 0;
+ M(c(), 0) = tmp0 > ForwardModel::Eps ? log(tmp0) : 0;
+ M(c(), 1) = tmp1 > ForwardModel::Eps ? log(tmp1) : 0;
if(FWD.GetProblem().GetNumL() > 0){
double tmp2 = Tmp(c(), 2);
- M(c(), 2) = tmp2> Eps ? log(tmp2) : 0;
+ M(c(), 2) = tmp2 > ForwardModel::Eps ? log(tmp2) : 0;
}
}
diff --git a/FWI/src/ForwardModels/ForwardModel.hpp b/FWI/src/ForwardModels/ForwardModel.hpp
index b0d690677cfc8d989308069d074a7102e3e6b5f4..c9f08bf21b17afda3b012361b863c5a6396056e9 100644
--- a/FWI/src/ForwardModels/ForwardModel.hpp
+++ b/FWI/src/ForwardModels/ForwardModel.hpp
@@ -144,6 +144,8 @@ protected:
MeasureComp);
}
public:
+ static constexpr double Eps = 1e-10;
+
std::unique_ptr<FWIObservationOperatorDGTS> Obs;
std::unique_ptr<FWITaper> Taper;
bool precon{false};
diff --git a/FWI/src/Helper/ABCHandler.cpp b/FWI/src/Helper/ABCHandler.cpp
index 1aa383c47dfbce4847b98d062e66af4c46cfa6a1..bef68e7e213a0c6c0804e1b8760705d6d7820c5f 100644
--- a/FWI/src/Helper/ABCHandler.cpp
+++ b/FWI/src/Helper/ABCHandler.cpp
@@ -27,11 +27,11 @@ void ABCHandler::findBndFlags(const Mesh &mesh,
for (int i = 0; i < dim; ++i) {
double d = abs(fp[i] - MeshBorders_[i * dim]);
double d2 = abs(fp[i] - MeshBorders_[i * dim + 1]);
- if (d < Eps) {
+ if (d < CosTaperConfig::Eps) {
bnds[i * dim] = bnd;
NrBndProcs[i * dim] = 1;
break;
- } else if (d2 < Eps) {
+ } else if (d2 < CosTaperConfig::Eps) {
bnds[i * dim + 1] = bnd;
NrBndProcs[i * dim + 1] = 1;
break;
@@ -121,7 +121,7 @@ void ABCHandler::initializeABCMesh(const Mesh &mesh,
std::vector<int> NrCells(dim, 0);
std::vector<double> PMLLength(2 * dim, 0.0);
auto meshWidth = mesh.MeshWidth();
- if (abs(meshWidth.first - meshWidth.second) > Eps)
+ if (abs(meshWidth.first - meshWidth.second) > CosTaperConfig::Eps)
THROW("Only for quadratic meshes")
dx = std::vector<double>(dim, meshWidth.first);
findCellNrAndABCLength(MeshBordersI, NrCells, NrBndProcs, PMLLength,
@@ -207,13 +207,15 @@ void ABCHandler::CreatePointsOutsideOfInnerGrid(
int ret = 0;
bool isOnLeftBnd = false;
for (; j < dim; ++j) {
- isOnLeftBnd = abs(fp[j] - MeshBorders[j * dim]) < Eps;
- bool isOnRightBnd = abs(fp[j] - MeshBorders[j * dim + 1]) < Eps;
+ isOnLeftBnd = abs(fp[j] - MeshBorders[j * dim]) < CosTaperConfig::Eps;
+ bool isOnRightBnd =
+ abs(fp[j] - MeshBorders[j * dim + 1]) < CosTaperConfig::Eps;
if (isOnLeftBnd || isOnRightBnd)
break;
}
if (abs(MeshBordersABC[j * dim + (isOnLeftBnd ? 0 : 1)] -
- MeshBorders[j * dim + (isOnLeftBnd ? 0 : 1)]) < Eps)
+ MeshBorders[j * dim + (isOnLeftBnd ? 0 : 1)]) <
+ CosTaperConfig::Eps)
continue;
double PMLLastBoundaryCellMid =
MeshBordersABC[j * dim + (isOnLeftBnd ? 0 : 1)] -
@@ -604,7 +606,8 @@ ABCMatVectorFWI::ABCMatVectorFWI(const ABCHandler &handler, const Vector &Unpadd
void ABCMatVectorFWI::applydistanceFactors(const Point &SendPoint, const Point &ReceivePoint, int i,
int j) {
- if (handler.GetDistFactor() < Eps) return;
+ if (handler.GetDistFactor() < CosTaperConfig::Eps)
+ return;
double distX = ReceivePoint[0] - SendPoint[0];
double distY = ReceivePoint[1] - SendPoint[1];
if (handler.GetConfig().getTaperStyle() == OLDTAPER){
@@ -770,22 +773,26 @@ ABCHandlerAAO::ABCHandlerAAO(const ABCConfigAAO &Conf_, const Vector &vecfull,in
};
for (cell c = vecfull.cells(); c!=vecfull.cells_end();++c){
position pos;
- if ((abs(c()[0] - dx[0]/2.0 - Conf.RectConf.LeftBnd)<Eps)){
+ if ((abs(c()[0] - dx[0] / 2.0 - Conf.RectConf.LeftBnd) <
+ CosTaperConfig::Eps)) {
if(c()[1] >Conf.RectConf.BottomBnd &&c()[1] < Conf.RectConf.TopBnd){
pos.left = true;
}
}
- if ((abs(c()[0] + dx[0]/2.0 - Conf.RectConf.RightBnd)<Eps)){
+ if ((abs(c()[0] + dx[0] / 2.0 - Conf.RectConf.RightBnd) <
+ CosTaperConfig::Eps)) {
if(c()[1] >Conf.RectConf.BottomBnd &&c()[1] < Conf.RectConf.TopBnd){
pos.right = true;
}
}
- if ((abs(c()[1] - dx[0]/2.0 - Conf.RectConf.BottomBnd)<Eps)){
+ if ((abs(c()[1] - dx[0] / 2.0 - Conf.RectConf.BottomBnd) <
+ CosTaperConfig::Eps)) {
if(c()[0] >Conf.RectConf.LeftBnd &&c()[0] < Conf.RectConf.RightBnd){
pos.bottom = true;
}
}
- if ((abs(c()[1] + dx[0]/2.0 - Conf.RectConf.TopBnd)<Eps)){
+ if ((abs(c()[1] + dx[0] / 2.0 - Conf.RectConf.TopBnd) <
+ CosTaperConfig::Eps)) {
if(c()[0] >Conf.RectConf.LeftBnd &&c()[0] < Conf.RectConf.RightBnd){
pos.top = true;
}
diff --git a/FWI/src/Helper/ABCHandler.hpp b/FWI/src/Helper/ABCHandler.hpp
index c24341a35138b644c21cd020b0b42d1a5b30c5af..925d17d9333ce2b651e41d9d1b2a02cbdf50cdcd 100644
--- a/FWI/src/Helper/ABCHandler.hpp
+++ b/FWI/src/Helper/ABCHandler.hpp
@@ -18,6 +18,9 @@ public:
double Qb = 1.1;
double alpha1 = -1;
double alpha2 = -2;
+
+ static constexpr double Eps = 1e-10;
+
CosTaperConfig(double alpha1_, double alpha2_, double deltad_ = 0.3, double vQ_ = 0.85,
double vB_ = 0.2, double Qb_ = 1.1) : alpha1(alpha1_),alpha2(alpha2_),
deltad(deltad_), vQ(vQ_),vB(vB_),Qb(Qb_){
@@ -345,12 +348,12 @@ public:
double distX, double distY, std::vector<double> LayerSizes) {
double relDist = -infty;
double Dist = -infty;
- if (abs(distX)>Eps) {
+ if (abs(distX) > CosTaperConfig::Eps) {
int direc = distX < 0 ? LEFT : RIGHT;
relDist = abs(distX) / LayerSizes[direc];
Dist = LayerSizes[direc];
}
- if (abs(distY) > Eps) {
+ if (abs(distY) > CosTaperConfig::Eps) {
int direc = distY < 0 ? DOWN : UP;
if (relDist < abs(distY) / LayerSizes[direc]){
relDist = abs(distY) / LayerSizes[direc];
@@ -361,7 +364,8 @@ public:
}
void applydistanceFactors(Vector& Vec,
const Point &SendPoint, const Point &ReceivePoint, int i) {
- if (DistFactor < Eps) return;
+ if (DistFactor < CosTaperConfig::Eps)
+ return;
double distX = ReceivePoint[0] - SendPoint[0];
double distY = ReceivePoint[1] - SendPoint[1];
if (Conf.getTaperStyle() == OLDTAPER){
@@ -438,11 +442,11 @@ class ABCMatVectorFWI : public ABCVectorFWI {
double calculateRelativeDistance(double distX,double distY){
double relDist = -infty;
- if (abs(distX)>Eps) {
+ if (abs(distX) > CosTaperConfig::Eps) {
int direc = distX < 0 ? LEFT : RIGHT;
relDist = abs(distX) / handler.GetConfig().GetABCPaddingLength(direc);
}
- if (abs(distY) > Eps) {
+ if (abs(distY) > CosTaperConfig::Eps) {
int direc = distY < 0 ? DOWN : UP;
relDist =
std::max(relDist,
diff --git a/FWI/src/Helper/FWIMppHelpers.cpp b/FWI/src/Helper/FWIMppHelpers.cpp
index 140102fb733f58ae4bb3bb9d7ab4e1589412a460..d0a75b267c20e92aa7b8d35badaaaf7a0933c605 100644
--- a/FWI/src/Helper/FWIMppHelpers.cpp
+++ b/FWI/src/Helper/FWIMppHelpers.cpp
@@ -45,6 +45,7 @@ double zKTauPositive(int numL, double alpha, double rhoK, double vK, double taup
std::vector<Point> PointsOnLine(std::vector<Point>& waypoints,
std::vector<int>& nrOfPointsOnLine) {
+ static constexpr double Eps = 1e-10;
unsigned int waypoint_count = waypoints.size();
std::vector<Point> PointsOnLine;
diff --git a/FWI/src/Helper/NormCalculation.hpp b/FWI/src/Helper/NormCalculation.hpp
index 4ae5979ef97bea741ad1387c0b345a5b2b4b540f..56d79938421e3d6700ff831d7f453740ef7e2fc5 100644
--- a/FWI/src/Helper/NormCalculation.hpp
+++ b/FWI/src/Helper/NormCalculation.hpp
@@ -12,8 +12,10 @@
#include "FWIHelpers.hpp"
class NormCalculation {
- const std::vector<unsigned int> &reconlist;
- const std::unordered_map<Point, DGAcousticElement *> &eleCache;
+ static constexpr double Eps = 1e-10;
+
+ const std::vector<unsigned int> &reconlist;
+ const std::unordered_map<Point, DGAcousticElement *> &eleCache;
public:
NormCalculation(const std::vector<unsigned int> &reconlist_,
const std::unordered_map<Point, DGAcousticElement *> &eleCache_
diff --git a/FWI/src/Observation/FWIObservationOperatorSpaceDGTS.hpp b/FWI/src/Observation/FWIObservationOperatorSpaceDGTS.hpp
index fc201e15c415b4212345822972c27da23a69e7a0..a65c5da192a97015d771d2596abdcf4106972932 100644
--- a/FWI/src/Observation/FWIObservationOperatorSpaceDGTS.hpp
+++ b/FWI/src/Observation/FWIObservationOperatorSpaceDGTS.hpp
@@ -75,6 +75,8 @@ class FWIObservationOperatorSpacePoint : public FWIObservationOperatorSpaceDGTS
void initalizeShapeValues(const ObservationSpecification &spec, const Mesh &b);
void initializeParallelizationData(const ObservationSpecification &spec){
+ static constexpr double Eps = 1e-10;
+
ParaDat = std::make_unique<SeismogramParallelizationData>(spec);
for (const auto& re : SpaceReceiver2SpaceSupport){
ParaDat->LocalReceivers.emplace_back(re.first);
diff --git a/FWI/src/Problems/FWIProblemsDGTSVA.cpp b/FWI/src/Problems/FWIProblemsDGTSVA.cpp
index 0610c90d18a5129455fc0e5391562b2ae2af8f9f..8ef3865f4967609271ca79f536dfc899da946f42 100644
--- a/FWI/src/Problems/FWIProblemsDGTSVA.cpp
+++ b/FWI/src/Problems/FWIProblemsDGTSVA.cpp
@@ -273,6 +273,8 @@ class FWIVAcousticSinCos2d_RHS : public WaveProblemDGTSVA {
Point taus;
public:
FWIVAcousticSinCos2d_RHS(int _numL = 1) : WaveProblemDGTSVA(_numL, 2) {
+ static constexpr double Eps = 1e-10;
+
setParametrization("native");
Config::Get("taus", taus);
Config::Get("kappa_rel", kappa_rel);
diff --git a/FWI/src/VAcousticMatFromVec.hpp b/FWI/src/VAcousticMatFromVec.hpp
index f87e4c048cb5c1219df6c529b20af30d2652eb98..231dcdc504b4439cfdd8f645bc080bd469bbb97d 100644
--- a/FWI/src/VAcousticMatFromVec.hpp
+++ b/FWI/src/VAcousticMatFromVec.hpp
@@ -11,7 +11,8 @@
#include "Problems/FWIAcousticInversionProb.hpp"
class VAcousticMatFromVec : public DGVAcousticMaterialFromVector {
- std::unique_ptr<Vector> TauP = nullptr;
+ static constexpr double Eps = 1e-10;
+ std::unique_ptr<Vector> TauP = nullptr;
const FWIVAcousticInversionProb &prob;
void fillVector(const FWIVAcousticInversionProb &problem_fwd, const Vector &Mat);
ABCHandler * handle=nullptr;
diff --git a/FWI/tests/FWITestLinearizedAdjoint.cpp b/FWI/tests/FWITestLinearizedAdjoint.cpp
index 0497310eded6c31bcc1ea0e3d87e146f878e95e4..2aa05976ed81eb529e8678750b68ff1a46fcb629 100644
--- a/FWI/tests/FWITestLinearizedAdjoint.cpp
+++ b/FWI/tests/FWITestLinearizedAdjoint.cpp
@@ -127,6 +127,8 @@ public:
};
TEST_P(FWITestLinearizedAdjoint, FWITestLinearizedAdjoint) {
+ static constexpr double Eps = 1e-10;
+
Vector Rand(*Mat);
Rand.Clear();
auto srcLoc = prob->getPointSrcLoc(0);