/** * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu * * See file AUTHORS for a list of contributors. * * This software is distributed under the terms of the GNU General Public * Licence version 3 (GPL Version 3). See file LICENSE for a full version of * the license. */ #define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one // cpp file #include <catch2/catch.hpp> #include <array> #include <iomanip> #include <iostream> #include <corsika/process/ProcessSequence.h> using namespace corsika; using namespace corsika::units::si; using namespace corsika::process; using namespace std; static const int nData = 10; int globalCount = 0; class ContinuousProcess1 : public ContinuousProcess<ContinuousProcess1> { int fV = 0; public: ContinuousProcess1(const int v) : fV(v) {} void Init() { cout << "ContinuousProcess1::Init" << endl; assert(globalCount == fV); globalCount++; } template <typename D, typename T, typename S> inline EProcessReturn DoContinuous(D& d, T&, S&) const { cout << "ContinuousProcess1::DoContinuous" << endl; for (int i = 0; i < nData; ++i) d.p[i] += 0.933; return EProcessReturn::eOk; } }; class ContinuousProcess2 : public ContinuousProcess<ContinuousProcess2> { int fV = 0; public: ContinuousProcess2(const int v) : fV(v) {} void Init() { cout << "ContinuousProcess2::Init" << endl; assert(globalCount == fV); globalCount++; } template <typename D, typename T, typename S> inline EProcessReturn DoContinuous(D& d, T&, S&) const { cout << "ContinuousProcess2::DoContinuous" << endl; for (int i = 0; i < nData; ++i) d.p[i] += 0.933; return EProcessReturn::eOk; } }; class Process1 : public DiscreteProcess<Process1> { public: Process1(const int v) : fV(v) {} void Init() { cout << "Process1::Init" << endl; assert(globalCount == fV); globalCount++; } template <typename D, typename S> inline EProcessReturn DoDiscrete(D& d, S&) const { for (int i = 0; i < nData; ++i) d.p[i] += 1 + i; return EProcessReturn::eOk; } // private: int fV; }; class Process2 : public DiscreteProcess<Process2> { int fV = 0; public: Process2(const int v) : fV(v) {} void Init() { cout << "Process2::Init" << endl; assert(globalCount == fV); globalCount++; } template <typename Particle, typename Stack> inline EProcessReturn DoDiscrete(Particle&, Stack&) const { cout << "Process2::DoDiscrete" << endl; return EProcessReturn::eOk; } template <typename Particle, typename Track> inline double GetInteractionLength(Particle&, Track&) const { cout << "Process2::GetInteractionLength" << endl; return 3; } }; class Process3 : public DiscreteProcess<Process3> { int fV = 0; public: Process3(const int v) : fV(v) {} void Init() { cout << "Process3::Init" << endl; assert(globalCount == fV); globalCount++; } template <typename Particle, typename Stack> inline EProcessReturn DoDiscrete(Particle&, Stack&) const { cout << "Process3::DoDiscrete" << endl; return EProcessReturn::eOk; } template <typename Particle, typename Track> inline double GetInteractionLength(Particle&, Track&) const { cout << "Process3::GetInteractionLength" << endl; return 1.; } }; class Process4 : public BaseProcess<Process4> { int fV = 0; public: Process4(const int v) : fV(v) {} void Init() { cout << "Process4::Init" << endl; assert(globalCount == fV); globalCount++; } template <typename D, typename T, typename S> inline EProcessReturn DoContinuous(D& d, T&, S&) const { for (int i = 0; i < nData; ++i) { d.p[i] /= 1.2; } return EProcessReturn::eOk; } // inline double MinStepLength(D& d) { template <typename Particle, typename Stack> EProcessReturn DoDiscrete(Particle&, Stack&) const { return EProcessReturn::eOk; } }; struct DummyData { double p[nData] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; }; struct DummyStack {}; struct DummyTrajectory {}; TEST_CASE("Process Sequence", "[Process Sequence]") { SECTION("Check init order") { Process1 m1(0); Process2 m2(1); Process3 m3(2); Process4 m4(3); const auto sequence = m1 + m2 + m3 + m4; globalCount = 0; sequence.Init(); // REQUIRE_NOTHROW( (sequence.Init()) ); // const auto sequence_wrong = m3 + m2 + m1 + m4; // globalCount = 0; // sequence_wrong.Init(); // REQUIRE_THROWS(sequence_wrong.Init()); } SECTION("interaction length") { ContinuousProcess1 cp1(0); Process2 m2(1); Process3 m3(2); DummyStack s; DummyTrajectory t; const auto sequence2 = cp1 + m2 + m3; double tot = sequence2.GetTotalInteractionLength(s, t); double tot_inv = sequence2.GetTotalInverseInteractionLength(s, t); cout << "lambda_tot=" << tot << " lambda_tot_inv=" << tot_inv << endl; } SECTION("sectionTwo") { ContinuousProcess1 cp1(0); ContinuousProcess2 cp2(3); Process2 m2(1); Process3 m3(2); const auto sequence2 = cp1 + m2 + m3 + cp2; DummyData p; DummyStack s; DummyTrajectory t; cout << "-->init sequence2" << endl; globalCount = 0; sequence2.Init(); cout << "-->docont" << endl; sequence2.DoContinuous(p, t, s); cout << "-->dodisc" << endl; sequence2.DoDiscrete(p, s); cout << "-->done" << endl; const int nLoop = 5; cout << "Running loop with n=" << nLoop << endl; for (int i = 0; i < nLoop; ++i) { sequence2.DoContinuous(p, t, s); sequence2.DoDiscrete(p, s); } for (int i = 0; i < nData; i++) { cout << "data[" << i << "]=" << p.p[i] << endl; } cout << "done" << endl; } }