Merge branch 'neumannseries' into 'master'

Neumannseries

See merge request kit/mpp/fwi!150

FWI/src/Helper/GramianBlockInverse.cpp

@@ -5,7 +5,12 @@
 #include "GramianBlockInverse.hpp"
 
 #include <utility>
-
+int truncate_to_8_significant_digits(double value) {
+  if (value == 0.0) { return 0.0; }
+  double magnitude = std::floor(std::log10(std::abs(value))) + 1;
+  double factor = std::pow(10, 8 - magnitude);
+  return int(std::trunc(value * factor));
+}
 RMatrix getRMatrixFromCell(const Matrix &A, const Cell &C) {
     std::vector<Point> z = A.GetDoF().GetStoragePoints(C);
     std::vector<row> R(z.size());
@@ -69,10 +74,10 @@ GramianBlockInverse::GramianBlockInverse(const Vector& WaveVec_,
         const DGAcousticElement *elem = elemIt.second;
         const Cell &C = elem->GetCell();
         double a = elem->Area();
-        if (ScalarMInverse.find(a) == ScalarMInverse.end()) {
-            RMatrix Sm = getRMatrixFromCell(*scalarProductMatrix, C);
+        if (ScalarMInverse.find(truncate_to_8_significant_digits(a)) == ScalarMInverse.end()) {
+          RMatrix Sm = getRMatrixFromCell(*scalarProductMatrix, C);
             Sm.Invert();
-            ScalarMInverse.emplace(a, SparseRMatrix(Sm));
+          ScalarMInverse.emplace(truncate_to_8_significant_digits(a), SparseRMatrix(Sm));
         }
     }
 }
@@ -88,7 +93,7 @@ void GramianBlockInverse::ScalarProductInverse(Vector &b) const {
 void GramianBlockInverse::productOneCellScalarProductInverse(const Vector &u, const DGAcousticElement &elem,
                                                              Vector &b) const {
     double a = elem.Area();
-  const SparseRMatrix &SpSm = ScalarMInverse.at(a);
+  const SparseRMatrix &SpSm = ScalarMInverse.at(truncate_to_8_significant_digits(a));
   int Row = u.find_row(elem.GetCell()()).Id();
   productOneCellSparse(u, SpSm, Row, b);
 }

FWI/src/Helper/GramianBlockInverse.hpp

@@ -20,8 +20,9 @@ class GramianBlockInverse {
     void productOneCellSparse(const Vector &u, const SparseRMatrix &SpSm, int Row, Vector &b) const;
 
     std::shared_ptr<Matrix> scalarProductMatrix = nullptr;
-    std::unordered_map<double, SparseRMatrix> ScalarMInverse;
-public:
+    std::unordered_map<int, SparseRMatrix> ScalarMInverse;
+
+  public:
     GramianBlockInverse(const Vector& WaveVec,
                         const std::unordered_map<Point, DGAcousticElement *> &eleCache);
 

FWI/src/Observation/DGSeismogram.cpp

@@ -131,9 +131,9 @@ void DGSeismogram::fillSeismogramWithNormalGaussianNoise(
   unsigned int localSeed = seed == 0 ? rd() : seed;
   unsigned int rnd = localSeed+PPM->Proc()*1000;
   std::normal_distribution<double> distribution(0.0, 1.0);
-  for (const auto& re : SeismoData.GetParaDat().LocalReceivers){
+  for (const auto &re : noiseSeismo.SeismoData.GetParaDat().LocalReceivers) {
     std::mt19937 generator(rnd++);
-    for (auto &t : GetTimeSeries(re)){
+    for (auto &t : noiseSeismo.GetTimeSeries(re)) {
       t = distribution(generator);
     }
   }