Line data    Source code

       1             : #include <stdexcept>
       2             : 
       3             : #include "crpropa/Grid.h"
       4             : #include "crpropa/Units.h"
       5             : #include "crpropa/Common.h"
       6             : #include "crpropa/GridTools.h"
       7             : #include "crpropa/magneticField/turbulentField/TurbulentField.h"
       8             : #include "crpropa/magneticField/turbulentField/GridTurbulence.h"
       9             : #include "crpropa/magneticField/turbulentField/PlaneWaveTurbulence.h"
      10             : #include "crpropa/magneticField/turbulentField/SimpleGridTurbulence.h"
      11             : 
      12             : #include "gtest/gtest.h"
      13             : 
      14             : using namespace crpropa;
      15             : 
      16             : //check problems brought up in https://github.com/CRPropa/CRPropa3/issues/322
      17           1 : TEST(testTurbulenceSpectrum, constructor) {
      18             :         double sIndex = 5./3.;
      19             :         double qIndex = 4.;
      20             :         double bendOver = 1.;
      21             :         double lMin = 1.; 
      22             :         double lMax = 10.; 
      23             :         double brms = 1*muG;
      24           1 :         auto spectrum = TurbulenceSpectrum(brms, lMin, lMax);
      25           1 :         EXPECT_DOUBLE_EQ(spectrum.getBrms(), brms);
      26           1 :         EXPECT_DOUBLE_EQ(spectrum.getLmin(), lMin);
      27           1 :         EXPECT_DOUBLE_EQ(spectrum.getLmax(), lMax);
      28           1 :         EXPECT_DOUBLE_EQ(spectrum.getLbendover(), bendOver); //default
      29           1 :         EXPECT_DOUBLE_EQ(spectrum.getSindex(), sIndex); //default
      30           1 :         EXPECT_DOUBLE_EQ(spectrum.getQindex(), qIndex); //default
      31           1 : }
      32             : 
      33           1 : TEST(testTurbulenceSpectrum, correlationLength) {
      34             :         double lMin = 0.00001; // not used for Lc
      35             :         double lMax = 9999999; // not used for Lc
      36             :         double lBo = 100;
      37           1 :         auto spectrum = TurbulenceSpectrum(1*muG, lMin, lMax, lBo);
      38             :         auto Lc = spectrum.getCorrelationLength();
      39           1 :     EXPECT_NEAR(Lc, 0.498*lBo, 0.001*lBo);
      40           1 : }
      41             : 
      42             : #ifdef CRPROPA_HAVE_FFTW3F
      43             : 
      44           1 : TEST(testSimpleGridTurbulence, oldFunctionForCrrelationLength) { //TODO: remove in future
      45             :         double lMin = 1*kpc;
      46             :         double lMax = 1*Gpc;
      47             :         double alpha = -11/3.;
      48           1 :         auto Lc = turbulentCorrelationLength(lMin, lMax, alpha);
      49           1 :     EXPECT_NEAR(Lc, lMax/5, 1*Mpc);
      50           1 : }
      51             : 
      52           2 : TEST(testVectorFieldGrid, Turbulence_bmean_brms) {
      53             :         // Test for zero mean: <B> = 0
      54             :         size_t n = 64;
      55             :         double spacing = 10 * Mpc / n;
      56             :         double Brms = 1;
      57             :         double lMin = 2 * spacing;
      58             :         double lMax = 8 * spacing;
      59             : 
      60             :         auto spectrum = SimpleTurbulenceSpectrum(Brms, lMin, lMax);
      61             :         auto gp = GridProperties(Vector3d(0, 0, 0), n, spacing);
      62           1 :     auto tf = SimpleGridTurbulence(spectrum, gp);
      63           1 :         auto grid = tf.getGrid();
      64             : 
      65             :         double precision = 1e-7;
      66           1 :         Vector3f bMean = meanFieldVector(grid);
      67           1 :         EXPECT_NEAR(0, bMean.x, precision);
      68           1 :         EXPECT_NEAR(0, bMean.y, precision);
      69           1 :         EXPECT_NEAR(0, bMean.z, precision);
      70           2 :         EXPECT_NEAR(1, rmsFieldStrength(grid), precision);
      71           1 : }
      72             : 
      73           2 : TEST(testVectorFieldGrid, Turbulence_seed) {
      74             :         // Test if seeding produces 2 identical fields
      75             :         size_t n = 64;
      76             :         double spacing = 1 * Mpc;
      77             :         double Brms = 1;
      78             :         double lMin = 2 * spacing;
      79             :         double lMax = 8 * spacing;
      80             :         int seed = 753;
      81             :         
      82             :         auto spectrum = SimpleTurbulenceSpectrum(Brms, lMin, lMax);
      83             : 
      84             :         auto gp1 = GridProperties(Vector3d(0, 0, 0), n, spacing);
      85           1 :     auto tf1 = SimpleGridTurbulence(spectrum, gp1,  seed);
      86             : 
      87             :         auto gp2 = GridProperties(Vector3d(0, 0, 0), n, spacing);
      88           1 :     auto tf2 = SimpleGridTurbulence(spectrum, gp2, seed);
      89             : 
      90             :         Vector3d pos(22 * Mpc);
      91           1 :         EXPECT_FLOAT_EQ(tf1.getField(pos).x, tf2.getField(pos).x);
      92           1 : }
      93             : 
      94             : #ifndef CRPROPA_TESTS_SKIP_EXCEPTIONS
      95           2 : TEST(testVectorFieldGrid, turbulence_Exceptions) {
      96             :         // Test exceptions
      97             :         size_t n = 64;
      98             :         double spacing = 10 * Mpc / n;
      99             :         double brms = 1;
     100           1 :         ref_ptr<Grid3f> grid = new Grid3f(Vector3d(0, 0, 0), n, spacing);
     101             : 
     102             :         // should be fine
     103           2 :         EXPECT_NO_THROW(initTurbulence(grid, brms, 2 * spacing, 8 * spacing));
     104             :         // lMin too small
     105           2 :         EXPECT_THROW(initTurbulence(grid, brms, 1.5 * spacing, 8 * spacing),
     106             :                         std::runtime_error);
     107             :         // lMin > lMax
     108           2 :         EXPECT_THROW(initTurbulence(grid, brms, 8.1 * spacing, 8 * spacing),
     109             :                         std::runtime_error);
     110             :         // lMax too large
     111           2 :         EXPECT_THROW(initTurbulence(grid, brms, 2 * spacing, 65 * spacing),
     112             :                         std::runtime_error);
     113           1 : }
     114             : #endif
     115             : 
     116           1 : TEST(testGridTurbulence, Turbulence_seed) {
     117             :         // Test if seeding produces 2 identical fields
     118             :         size_t n = 64;
     119             :         double spacing = 1 * Mpc;
     120             :         double Brms = 1;
     121             :         double lMin = 2 * spacing;
     122             :         double lMax = 8 * spacing;
     123             :         double lBo = lMax/6;
     124             :         int seed = 137;
     125             :         
     126           1 :         auto spectrum = TurbulenceSpectrum(Brms, lMin, lMax, lBo);
     127             : 
     128             :         auto gp1 = GridProperties(Vector3d(0, 0, 0), n, spacing);
     129           1 :     auto tf1 = GridTurbulence(spectrum, gp1, seed);
     130             : 
     131             :         auto gp2 = GridProperties(Vector3d(0, 0, 0), n, spacing);
     132           1 :     auto tf2 = GridTurbulence(spectrum, gp2, seed);
     133             : 
     134             :         Vector3d pos(22 * Mpc);
     135           1 :         EXPECT_FLOAT_EQ(tf1.getField(pos).x, tf2.getField(pos).x);
     136           1 : }
     137             : #endif // CRPROPA_HAVE_FFTW3F
     138             : 
     139           1 : int main(int argc, char **argv) {
     140           1 :         ::testing::InitGoogleTest(&argc, argv);
     141             :         return RUN_ALL_TESTS();
     142             : }