Merge branch 'develop' into MultiThreadV2

Core/Multilayer/SpecularScalarStrategy.cpp

@@ -95,21 +95,17 @@ bool SpecularScalarStrategy::calculateUpFromLayer(std::vector<ScalarRTCoefficien
 
         coeff[i].t_r = transition(kz[i], kz[i + 1], sigma, slices[i].thickness(), coeff[i + 1].t_r);
 
-        // normalize the t-coefficient in each step of the computation
-        // this avoids an overflow in the backwards propagation
-        // the used factors are stored in order to correct the amplitudes
-        // in forward direction at the end of the computation
-
-        if (std::isinf(std::norm(coeff[i].t_r(0)))){
-//            throw std::runtime_error("Unexpected inf");
-            // catching this inf is not correctly covered by existing tests
-            // it is possible to pass all tests by moving the normalization before
-            // this check and hence corrupting all amplitudes
+        if (std::isinf(std::norm(coeff[i].t_r(0))) || std::isnan(std::norm(coeff[i].t_r(0)))){
             coeff[i].t_r(0) = 1.0;
             coeff[i].t_r(1) = 0.0;
 
             setZeroBelow(coeff, i);
         }
+
+        // normalize the t-coefficient in each step of the computation
+        // this avoids an overflow in the backwards propagation
+        // the used factors are stored in order to correct the amplitudes
+        // in forward direction at the end of the computation
         factors[i] = coeff[i].t_r(0);
         coeff[i].t_r = coeff[i].t_r / coeff[i].t_r(0);
     }

Core/StandardSamples/MultiLayerWithNCRoughnessBuilder.cpp

+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Core/StandardSamples/MultiLayerWithNCRoughnessBuilder.cpp
+//! @brief     Implement class MultiLayerWithNCRoughnessBuilder.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************** //
+
+#include "MultiLayerWithNCRoughnessBuilder.h"
+#include "MaterialFactoryFuncs.h"
+#include "Layer.h"
+#include "LayerRoughness.h"
+#include "MultiLayer.h"
+#include "RealParameter.h"
+#include "Units.h"
+#include "BornAgainNamespace.h"
+
+MultiLayerWithNCRoughnessBuilder::MultiLayerWithNCRoughnessBuilder()
+{}
+
+MultiLayer* MultiLayerWithNCRoughnessBuilder::buildSample() const
+{
+    auto multi_layer = MultiLayerWithRoughnessBuilder::buildSample();
+    multi_layer->setRoughnessModel(RoughnessModel::NEVOT_CROCE);
+    return multi_layer;
+}

Core/StandardSamples/MultiLayerWithNCRoughnessBuilder.h

+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Core/StandardSamples/MultiLayerWithNCRoughnessBuilder.h
+//! @brief     Defines class MultiLayerWithNCRoughnessBuilder.
+//!
+//! @homepage  http://www.bornagainproject.org
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2018
+//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
+//
+// ************************************************************************** //
+
+#ifndef MULTILAYERWITHNCROUGHNESSBUILDER_H
+#define MULTILAYERWITHNCROUGHNESSBUILDER_H
+
+#include "IMultiLayerBuilder.h"
+#include "MultiLayerWithRoughnessBuilder.h"
+
+class ISample;
+
+//! Builds sample: layers with correlated roughness.
+//! @ingroup standard_samples
+
+class BA_CORE_API_ MultiLayerWithNCRoughnessBuilder : public MultiLayerWithRoughnessBuilder
+{
+public:
+    MultiLayerWithNCRoughnessBuilder();
+    MultiLayer* buildSample() const override;
+
+private:
+};
+
+#endif // MULTILAYERWITHNCROUGHNESSBUILDER_H

Core/StandardSamples/MultiLayerWithRoughnessBuilder.h

@@ -26,7 +26,7 @@ class BA_CORE_API_ MultiLayerWithRoughnessBuilder : public IMultiLayerBuilder
 {
 public:
     MultiLayerWithRoughnessBuilder();
-    MultiLayer* buildSample() const;
+    virtual MultiLayer* buildSample() const;
 
 private:
     double m_thicknessA;

Core/StandardSamples/SampleBuilderFactory.cpp

@@ -27,6 +27,7 @@
 #include "MagneticParticlesBuilder.h"
 #include "MesoCrystalBuilder.h"
 #include "MultiLayerWithRoughnessBuilder.h"
+#include "MultiLayerWithNCRoughnessBuilder.h"
 #include "MultipleLayoutBuilder.h"
 #include "ParaCrystalBuilder.h"
 #include "ParticleCompositionBuilder.h"
@@ -152,6 +153,9 @@ SampleBuilderFactory::SampleBuilderFactory()
     registerItem("MultiLayerWithRoughnessBuilder", create_new<MultiLayerWithRoughnessBuilder>,
                  "Layer with correlated roughness");
 
+    registerItem("MultiLayerWithNCRoughnessBuilder", create_new<MultiLayerWithNCRoughnessBuilder>,
+                 "Layer with correlated roughness");
+    
     registerItem("TwoLayerRoughnessBuilder", create_new<TwoLayerRoughnessBuilder>,
                  "Two layers with rough interface");
 

Tests/Functional/Core/CoreStandardTest/CMakeLists.txt

@@ -74,6 +74,7 @@ set(test_cases
     HomogeneousTiNiSample
     HomogeneousTiNiSampleWithAbsorption
     RoughnessInSpecular
+    NCRoughnessInSpecular
     GaussianBeamFootprint
     SquareBeamFootprint
     SpecularDivergentBeam

Tests/Functional/TestMachinery/StandardTestCatalogue.cpp

@@ -465,6 +465,12 @@ StandardTestCatalogue::StandardTestCatalogue()
         "MultiLayerWithRoughnessBuilder",
         2e-9);
 
+    add("NCRoughnessInSpecular",
+        "Specular simulation with rough sample (Nevot-Croce)",
+        "BasicSpecular",
+        "MultiLayerWithNCRoughnessBuilder",
+        2e-9);
+    
     add("GaussianBeamFootprint",
         "Similar to HomogeneousTiNiSample, but with finite-sized gaussian beam",
         "SpecularWithGaussianBeam",

Tests/UnitTests/Core/Sample/RTTest.cpp

@@ -46,7 +46,7 @@ protected:
     const Material stone = HomogeneousMaterial("substrate material", 1e-6, 1e-7);
     const Layer topLayer{air, 0};
     const Layer substrate{stone, 0};
-    const kvector_t k{1, 0, 1e-3};
+    const kvector_t k{1, 0, -1e-3};
     MultiLayer sample1, sample2;
     std::vector<ScalarRTCoefficients> coeffs1, coeffs2;
 };
@@ -110,8 +110,50 @@ TEST_F(RTTest, SplitBilayers)
     coeffs1 = getCoeffs( std::make_unique<SpecularScalarTanhStrategy>()->Execute(sample_1.slices(), k) );
     coeffs2 = getCoeffs( std::make_unique<SpecularScalarTanhStrategy>()->Execute(sample_2.slices(), k) );
 
-    // printCoeffs( coeffs1 );
-    // printCoeffs( coeffs2 );
+//     printCoeffs( coeffs1 );
+//     printCoeffs( coeffs2 );
+
+    compareCoeffs(coeffs1[0], coeffs2[0]);
+    compareCoeffs(coeffs1[1], coeffs2[1]);
+
+    // Amplitudes at bottom must be strictly zero.
+    // The new algorithm handles this without an overflow
+    EXPECT_EQ(complex_t(), coeffs1[coeffs1.size() - 2].t_r(0));
+    EXPECT_EQ(complex_t(), coeffs1[coeffs1.size() - 2].t_r(1));
+    EXPECT_EQ(complex_t(), coeffs2[coeffs2.size() - 2].t_r(0));
+    EXPECT_EQ(complex_t(), coeffs2[coeffs2.size() - 2].t_r(1));
+}
+
+TEST_F(RTTest, Overflow)
+{
+    // Text extra thick layers to also provoke an overflow in the new algorithm
+    const int n = 5;
+
+    sample1.addLayer(topLayer);
+    for (size_t i = 0; i < n; ++i) {
+        sample1.addLayer(Layer(amat, 1000));
+        sample1.addLayer(Layer(bmat, 200000));
+    }
+    sample1.addLayer(substrate);
+
+    sample2.addLayer(topLayer);
+    for (size_t i = 0; i < n; ++i) {
+        sample2.addLayer(Layer(amat, 1000));
+        sample2.addLayer(Layer(bmat, 100000));
+        sample2.addLayer(Layer(amat, 0));
+        sample2.addLayer(Layer(bmat, 100000));
+    }
+    sample2.addLayer(substrate);
+
+    SimulationOptions options;
+    ProcessedSample sample_1(sample1, options);
+    ProcessedSample sample_2(sample2, options);
+
+    coeffs1 = getCoeffs( std::make_unique<SpecularScalarTanhStrategy>()->Execute(sample_1.slices(), k) );
+    coeffs2 = getCoeffs( std::make_unique<SpecularScalarTanhStrategy>()->Execute(sample_2.slices(), k) );
+
+//     printCoeffs( coeffs1 );
+//     printCoeffs( coeffs2 );
 
     compareCoeffs(coeffs1[0], coeffs2[0]);
     compareCoeffs(coeffs1[1], coeffs2[1]);