From 159919b19f88fdf9f144ed3730c4b9e951f7a9cd Mon Sep 17 00:00:00 2001 From: ralfulrich <ralf.ulrich@kit.edu> Date: Tue, 10 Nov 2020 09:52:26 +0100 Subject: [PATCH] fixed testCascade --- .gitlab-ci.yml | 17 ++--- Environment/CMakeLists.txt | 1 + Environment/Environment.h | 1 - Environment/VolumeTreeNode.h | 11 ++-- Framework/Cascade/Cascade.h | 13 ++-- Framework/Cascade/testCascade.cc | 64 ++++++++++++------- Framework/Cascade/testCascade.h | 5 +- Framework/ProcessSequence/ProcessSequence.h | 10 +-- Processes/Tracking/Intersect.hpp | 53 +++++++-------- Processes/TrackingLeapFrogCurved/Tracking.h | 29 +++++++-- Processes/TrackingLeapFrogStraight/Tracking.h | 4 +- 11 files changed, 113 insertions(+), 95 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 464ce2402..9ebe0555b 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -170,7 +170,7 @@ build-clang-8: script: - cd build - set -o pipefail - - ctest -j4 -VV | gzip -v -9 > test.log.gz + - ctest -j4 rules: - if: $CI_MERGE_REQUEST_ID when: manual @@ -185,8 +185,6 @@ build-clang-8: reports: junit: - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml - paths: - - ${CI_PROJECT_DIR}/build/test.log.gz cache: paths: - ${CI_PROJECT_DIR}/build/ @@ -227,7 +225,7 @@ test-clang-8: - cd build - cmake --build . -- -j4 - set -o pipefail - - ctest -j4 -VV | gzip -v -9 > test.log.gz + - ctest -j4 rules: - if: '$CI_MERGE_REQUEST_ID && $CI_MERGE_REQUEST_TITLE =~ /^Draft:/' allow_failure: false @@ -241,8 +239,6 @@ test-clang-8: reports: junit: - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml - paths: - - ${CI_PROJECT_DIR}/build/test.log.gz cache: paths: - ${CI_PROJECT_DIR}/build/ @@ -338,7 +334,7 @@ example-clang-8: - cd build - cmake --build . -- -j4 - set -o pipefail - - ctest -j4 -VV | gzip -v -9 > test.log.gz + - ctest -j4 - make -j4 run_examples | gzip -v -9 > examples.log.gz rules: - if: '$CI_MERGE_REQUEST_ID && $CI_MERGE_REQUEST_TITLE =~ /^Draft:/' @@ -357,7 +353,6 @@ example-clang-8: - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml paths: - ${CI_PROJECT_DIR}/build/examples.log.gz - - ${CI_PROJECT_DIR}/build/test.log.gz cache: paths: - ${CI_PROJECT_DIR}/build/ @@ -450,7 +445,7 @@ install-clang-8: - cmake .. -DCMAKE_BUILD_TYPE=Release - cmake --build . -- -j4 - set -o pipefail - - ctest -j4 -VV | gzip -v -9 > test.log.gz + - ctest -j4 - make -j4 run_examples | gzip -v -9 > examples.log.gz rules: - if: '$CI_MERGE_REQUEST_LABELS =~ /Ready for code review/' # run on merge requests, if label 'Ready for code review' is set @@ -475,7 +470,6 @@ install-clang-8: junit: - ${CI_PROJECT_DIR}/build/test_outputs/junit*.xml paths: - - ${CI_PROJECT_DIR}/build/test.log.gz - ${CI_PROJECT_DIR}/build/examples.log.gz # release for gcc @@ -513,7 +507,7 @@ coverage: - cd build - cmake .. -DCMAKE_BUILD_TYPE=Coverage - cmake --build . -- -j4 - - ctest -j4 -VV | gzip -v -9 > test.log.gz + - ctest -j4 - cmake --build . --target coverage - tar czf coverage-report.tar.gz coverage-report coverage: '/^.*functions\.+:\s(.*\%)\s/' @@ -530,7 +524,6 @@ coverage: when: always expire_in: 1 year paths: - - ${CI_PROJECT_DIR}/build/test.log.gz - ${CI_PROJECT_DIR}/build/coverage-report.tar.gz cache: paths: diff --git a/Environment/CMakeLists.txt b/Environment/CMakeLists.txt index 46ecfb3cb..3231ab4a8 100644 --- a/Environment/CMakeLists.txt +++ b/Environment/CMakeLists.txt @@ -19,6 +19,7 @@ set ( set ( ENVIRONMENT_HEADERS VolumeTreeNode.h + IEmpty.hpp IMediumModel.h NuclearComposition.h HomogeneousMedium.h diff --git a/Environment/Environment.h b/Environment/Environment.h index 774b74927..1f807de45 100644 --- a/Environment/Environment.h +++ b/Environment/Environment.h @@ -8,7 +8,6 @@ #pragma once -#include <corsika/environment/IMediumModel.h> #include <corsika/environment/VolumeTreeNode.h> #include <corsika/geometry/Point.h> #include <corsika/geometry/RootCoordinateSystem.h> diff --git a/Environment/VolumeTreeNode.h b/Environment/VolumeTreeNode.h index 8cad83555..ad3306810 100644 --- a/Environment/VolumeTreeNode.h +++ b/Environment/VolumeTreeNode.h @@ -8,16 +8,14 @@ #pragma once -#include <corsika/environment/IMediumModel.h> +#include <corsika/environment/IEmpty.hpp> #include <corsika/geometry/Volume.h> #include <memory> #include <vector> namespace corsika::environment { - class Empty {}; //<! intended for usage as default template argument - - template <typename TModelProperties = Empty> + template <typename TModelProperties = IEmpty> class VolumeTreeNode { public: using IModelProperties = TModelProperties; @@ -116,14 +114,15 @@ namespace corsika::environment { void SetModelProperties(IMPSharedPtr ptr) { fModelProperties = ptr; } + /* template <class MediumType, typename... Args> static auto CreateMedium(Args&&... args) { - static_assert(std::is_base_of_v<IMediumModel, MediumType>, + static_assert(std::is_base_of_v<, MediumType>, "unusable type provided, needs to be derived from \"IMediumModel\""); return std::make_shared<MediumType>(std::forward<Args>(args)...); } - + */ private: std::vector<VTNUPtr> fChildNodes; std::vector<VolumeTreeNode<IModelProperties> const*> fExcludedNodes; diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h index 3e81ef973..6dc602de2 100644 --- a/Framework/Cascade/Cascade.h +++ b/Framework/Cascade/Cascade.h @@ -19,8 +19,6 @@ #include <corsika/stack/history/EventType.hpp> #include <corsika/stack/history/HistorySecondaryProducer.hpp> -#include <corsika/setup/SetupTrajectory.h> - /* see Issue 161, we need to include SetupStack only because we need to globally define StackView. This is clearly not nice and should be changed, when possible. It might be that StackView needs to be @@ -174,8 +172,9 @@ namespace corsika::cascade { // assert that particle stays outside void Universe if it has no // model properties set - assert(currentLogicalNode != &*environment_.GetUniverse() || - environment_.GetUniverse()->HasModelProperties()); + assert((currentLogicalNode != &*environment_.GetUniverse() || + environment_.GetUniverse()->HasModelProperties()) && + "FATAL: The environment model has no valid properties set!"); // determine combined total inverse decay time InverseTimeType const total_inv_lifetime = @@ -237,9 +236,9 @@ namespace corsika::cascade { } C8LOG_DEBUG("sth. happening before geometric limit ? {}", - ((min_distance < geomMaxLength) ? "yes" : "no")); + ((min_distance <= geomMaxLength) ? "yes" : "no")); - if (min_distance < geomMaxLength) { // interaction to happen within geometric limit + if (min_distance <= geomMaxLength) { // interaction to happen within geometric limit // check whether decay or interaction limits this step the // outcome of decay or interaction MAY be a) new particles in @@ -333,7 +332,7 @@ namespace corsika::cascade { const auto sample_process = uniDist(rng_); auto const returnCode = process_sequence_.SelectInteraction(view, sample_process); if (returnCode != process::EProcessReturn::eInteracted) { - C8LOG_WARN("Particle did not interace!"); + C8LOG_WARN("Particle did not interact!"); } SetEventType(view, history::EventType::Interaction); return returnCode; diff --git a/Framework/Cascade/testCascade.cc b/Framework/Cascade/testCascade.cc index 1018bb002..2ca4dd7f8 100644 --- a/Framework/Cascade/testCascade.cc +++ b/Framework/Cascade/testCascade.cc @@ -13,8 +13,6 @@ #include <corsika/process/ProcessSequence.h> #include <corsika/process/NullModel.h> #include <corsika/process/stack_inspector/StackInspector.h> -#include <corsika/process/tracking_line/Tracking.h> -//#include <corsika/process/tracking_bfield/Tracking.h> #include <corsika/particles/ParticleProperties.h> @@ -36,44 +34,61 @@ using namespace corsika::geometry; #include <limits> using namespace std; -/* - The dummy env must support GetMagneticField(), and a density model +/** + * testCascade implements an e.m. Heitler model with energy splitting + * and a critical energy. + * + * It resembles one of the most simple cascades you can simulate with CORSIKA8. + **/ +/* + The dummy env (here) doesn't need to have any propoerties */ auto MakeDummyEnv() { TestEnvironmentType env; // dummy environment auto& universe = *(env.GetUniverse()); - const geometry::CoordinateSystem& cs = env.GetCoordinateSystem(); - auto theMedium = TestEnvironmentType::CreateNode<Sphere>( + auto world = TestEnvironmentType::CreateNode<Sphere>( Point{env.GetCoordinateSystem(), 0_m, 0_m, 0_m}, 1_m * std::numeric_limits<double>::infinity()); - using MyHomogeneousModel = environment::UniformMagneticField< - environment::HomogeneousMedium<TestEnvironmentInterface>>; + using MyEmptyModel = environment::Empty<environment::IEmpty>; + world->SetModelProperties<MyEmptyModel>(); - theMedium->SetModelProperties<MyHomogeneousModel>( - geometry::Vector(cs, 0_T, 0_T, 0_T), 1_g / (1_cm * 1_cm * 1_cm), - environment::NuclearComposition( - std::vector<particles::Code>{particles::Code::Proton}, std::vector<float>{1.})); - - universe.AddChild(std::move(theMedium)); + universe.AddChild(std::move(world)); return env; } +/** + * \class DummyTracking + * + * For the Heitler model we don't need particle transport. + **/ +class DummyTracking { + +public: + template <typename TParticle> + auto GetTrack(TParticle const& particle) { + using namespace corsika::units::si; + using namespace corsika::geometry; + geometry::Vector<SpeedType::dimension_type> const initialVelocity = + particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c; + return std::make_tuple( + geometry::LineTrajectory(geometry::Line(particle.GetPosition(), initialVelocity), + 0_s), // trajectory + particle.GetNode()); // next volume node + } +}; + class ProcessSplit : public process::InteractionProcess<ProcessSplit> { int fCalls = 0; - GrammageType fX0; public: - ProcessSplit(GrammageType const X0) - : fX0(X0) {} - template <typename Particle> corsika::units::si::GrammageType GetInteractionLength(Particle const&) const { - return fX0; + return 0_g / square(1_cm); } template <typename TSecondaryView> @@ -126,6 +141,8 @@ public: TEST_CASE("Cascade", "[Cascade]") { + logging::SetLevel(logging::level::trace); + HEPEnergyType E0 = 100_GeV; random::RNGManager& rmng = random::RNGManager::GetInstance(); @@ -133,14 +150,12 @@ TEST_CASE("Cascade", "[Cascade]") { auto env = MakeDummyEnv(); auto const& rootCS = env.GetCoordinateSystem(); - tracking_line::Tracking tracking; stack_inspector::StackInspector<TestCascadeStack> stackInspect(1, true, E0); process::NullModel nullModel; - const GrammageType X0 = 20_g / square(1_cm); const HEPEnergyType Ecrit = 85_MeV; - ProcessSplit split(X0); + ProcessSplit split; ProcessCut cut(Ecrit); auto sequence = process::sequence(nullModel, stackInspect, split, cut); TestCascadeStack stack; @@ -153,7 +168,8 @@ TEST_CASE("Cascade", "[Cascade]") { particles::GetMass(particles::Code::Electron)))}), Point(rootCS, {0_m, 0_m, 10_km}), 0_ns)); - cascade::Cascade<tracking_line::Tracking, decltype(sequence), TestCascadeStack, + DummyTracking tracking; + cascade::Cascade<DummyTracking, decltype(sequence), TestCascadeStack, TestCascadeStackView> EAS(env, tracking, sequence, stack); @@ -161,7 +177,7 @@ TEST_CASE("Cascade", "[Cascade]") { EAS.Run(); CHECK(cut.GetCount() == 2048); - CHECK(cut.GetCalls() == 2047); + CHECK(cut.GetCalls() == 2047); // final particle is still on stack and not yet deleted CHECK(split.GetCalls() == 2047); } diff --git a/Framework/Cascade/testCascade.h b/Framework/Cascade/testCascade.h index 848c597ee..e05744143 100644 --- a/Framework/Cascade/testCascade.h +++ b/Framework/Cascade/testCascade.h @@ -9,15 +9,14 @@ #pragma once #include <corsika/environment/Environment.h> -#include <corsika/environment/IMagneticFieldModel.h> +#include <corsika/environment/IEmpty.hpp> #include <corsika/stack/CombinedStack.h> #include <corsika/stack/SecondaryView.h> #include <corsika/stack/node/GeometryNodeStackExtension.h> #include <corsika/stack/nuclear_extension/NuclearStackExtension.h> -using TestEnvironmentInterface = - corsika::environment::IMagneticFieldModel<corsika::environment::IMediumModel>; +using TestEnvironmentInterface = corsika::environment::IEmpty; using TestEnvironmentType = corsika::environment::Environment<TestEnvironmentInterface>; template <typename T> diff --git a/Framework/ProcessSequence/ProcessSequence.h b/Framework/ProcessSequence/ProcessSequence.h index 32d302594..02a3d062e 100644 --- a/Framework/ProcessSequence/ProcessSequence.h +++ b/Framework/ProcessSequence/ProcessSequence.h @@ -215,7 +215,7 @@ namespace corsika::process { const auto particle = view.parent(); lambda_inv_sum += A_.GetInverseInteractionLength(particle); // check if we should execute THIS process and then EXIT - if (lambda_inv_select < lambda_inv_sum) { + if (lambda_inv_select <= lambda_inv_sum) { A_.DoInteraction(view); return EProcessReturn::eInteracted; } @@ -229,7 +229,7 @@ namespace corsika::process { // if this is not a ContinuousProcess --> evaluate probability lambda_inv_sum += B_.GetInverseInteractionLength(view.parent()); // check if we should execute THIS process and then EXIT - if (lambda_inv_select < lambda_inv_sum) { + if (lambda_inv_select <= lambda_inv_sum) { B_.DoInteraction(view); return EProcessReturn::eInteracted; } @@ -279,8 +279,8 @@ namespace corsika::process { // if this is not a ContinuousProcess --> evaluate probability decay_inv_sum += A_.GetInverseLifetime(view.parent()); // check if we should execute THIS process and then EXIT - if (decay_inv_select < decay_inv_sum) { // more pedagogical: rndm_select < - // decay_inv_sum / decay_inv_tot + if (decay_inv_select <= decay_inv_sum) { // more pedagogical: rndm_select < + // decay_inv_sum / decay_inv_tot A_.DoDecay(view); return EProcessReturn::eDecayed; } @@ -294,7 +294,7 @@ namespace corsika::process { // if this is not a ContinuousProcess --> evaluate probability decay_inv_sum += B_.GetInverseLifetime(view.parent()); // check if we should execute THIS process and then EXIT - if (decay_inv_select < decay_inv_sum) { + if (decay_inv_select <= decay_inv_sum) { B_.DoDecay(view); return EProcessReturn::eDecayed; } diff --git a/Processes/Tracking/Intersect.hpp b/Processes/Tracking/Intersect.hpp index 277b66f02..277c5f633 100644 --- a/Processes/Tracking/Intersect.hpp +++ b/Processes/Tracking/Intersect.hpp @@ -14,12 +14,25 @@ #include <corsika/logging/Logging.h> #include <corsika/geometry/Intersections.hpp> +#include <limits> + namespace corsika::process::tracking { /** * \class Intersect * + * This is a CRTP class to provide a generic volume-tree + * intersection for the purpose of tracking. + * + * It return the closest distance in time to the next geometric + * intersection, as well as a pointer to the corresponding new + * volume. * + * User may provide an optional global step-length limit as + * parameter. This may be needd for (simpler) algorithms in magnetic + * fields, where tracking errors grow linearly with step-length. + * Obviously, in the case of the step-length limit, the returend + * "next" volume is just the current one. * **/ @@ -28,45 +41,23 @@ namespace corsika::process::tracking { protected: template <typename TParticle> - auto nextIntersect(const TParticle& particle) const { + auto nextIntersect( + const TParticle& particle, + corsika::units::si::TimeType step_limit = + std::numeric_limits<corsika::units::si::TimeType::value_type>::infinity() * + corsika::units::si::second) const { using namespace corsika::units::si; using namespace corsika::geometry; - const Point& initialPosition = particle.GetPosition(); - typedef typename std::remove_reference<decltype(*particle.GetNode())>::type node_type; - node_type& volumeNode = *particle.GetNode(); + node_type& volumeNode = + *particle.GetNode(); // current "logical" node, from previous tracking step C8LOG_DEBUG("volumeNode={}, numericallyInside={} ", fmt::ptr(&volumeNode), - volumeNode.GetVolume().Contains(initialPosition)); - - auto const velocity = - particle.GetMomentum() / particle.GetEnergy() * corsika::units::constants::c; - - // for the event of magnetic fields and curved trajectories, we need to limit - // maximum step-length since we need to follow curved - // trajectories segment-wise -- at least if we don't employ concepts as "Helix - // Trajectories" or similar - const auto& magneticfield = - volumeNode.GetModelProperties().GetMagneticField(initialPosition); - const auto magnitudeB = magneticfield.norm(); - const int chargeNumber = particle.GetChargeNumber(); - auto const momentumVerticalMag = - particle.GetMomentum() - - particle.GetMomentum().parallelProjectionOnto(magneticfield); - LengthType const gyroradius = - (chargeNumber == 0 || magnitudeB == 0_T - ? std::numeric_limits<TimeType::value_type>::infinity() * 1_m - : momentumVerticalMag.norm() * 1_V / - (corsika::units::constants::c * abs(chargeNumber) * magnitudeB * - 1_eV)); - const double maxRadians = 0.01; - const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius; - C8LOG_DEBUG("gyroradius {}, steplimit: {} m = {} s", gyroradius, steplimit, - steplimit / velocity.norm()); + volumeNode.GetVolume().Contains(particle.GetPosition())); // start values: - TimeType minTime = steplimit / velocity.norm(); + TimeType minTime = step_limit; node_type* minNode = &volumeNode; // determine the first geometric collision with any other Volume boundary diff --git a/Processes/TrackingLeapFrogCurved/Tracking.h b/Processes/TrackingLeapFrogCurved/Tracking.h index ee0ba6e7f..4ffa43828 100644 --- a/Processes/TrackingLeapFrogCurved/Tracking.h +++ b/Processes/TrackingLeapFrogCurved/Tracking.h @@ -107,13 +107,34 @@ namespace corsika::process { typename std::remove_reference<decltype(*particle.GetNode())>::type node_type; node_type& volumeNode = *particle.GetNode(); + // for the event of magnetic fields and curved trajectories, we need to limit + // maximum step-length since we need to follow curved + // trajectories segment-wise -- at least if we don't employ concepts as "Helix + // Trajectories" or similar + const auto& magneticfield = + volumeNode.GetModelProperties().GetMagneticField(position); + const auto magnitudeB = magneticfield.norm(); + const int chargeNumber = particle.GetChargeNumber(); + auto const momentumVerticalMag = + particle.GetMomentum() - + particle.GetMomentum().parallelProjectionOnto(magneticfield); + LengthType const gyroradius = + (chargeNumber == 0 || magnitudeB == 0_T + ? std::numeric_limits<TimeType::value_type>::infinity() * 1_m + : momentumVerticalMag.norm() * 1_V / + (corsika::units::constants::c * abs(chargeNumber) * magnitudeB * + 1_eV)); + const double maxRadians = 0.01; + const LengthType steplimit = 2 * cos(maxRadians) * sin(maxRadians) * gyroradius; + const TimeType steplimit_time = steplimit / initialVelocity.norm(); + C8LOG_DEBUG("gyroradius {}, steplimit: {} m = {} s", gyroradius, steplimit, + steplimit_time); + // traverse the environment volume tree and find next // intersection - auto [minTime, minNode] = tracking::Intersect<Tracking>::nextIntersect(particle); + auto [minTime, minNode] = + tracking::Intersect<Tracking>::nextIntersect(particle, steplimit_time); - const int chargeNumber = particle.GetChargeNumber(); - const auto& magneticfield = - volumeNode.GetModelProperties().GetMagneticField(position); const auto k = chargeNumber * corsika::units::constants::cSquared * 1_eV / (particle.GetEnergy() * 1_V); return std::make_tuple( diff --git a/Processes/TrackingLeapFrogStraight/Tracking.h b/Processes/TrackingLeapFrogStraight/Tracking.h index fd15b89c8..d9e7c57de 100644 --- a/Processes/TrackingLeapFrogStraight/Tracking.h +++ b/Processes/TrackingLeapFrogStraight/Tracking.h @@ -39,14 +39,14 @@ namespace corsika::process { * tracking_line::Tracking and adds a (two-step) Leap-Frog * algorithms with two halve-steps and magnetic deflection. * - * The two halve steps are implemented as two + * The two halve steps are implemented as two straight explicit * `tracking_line::Tracking`s and all geometry intersections are, * thus, based on those two straight line elements. * * As a precaution for numerical instability, the steplength is * limited to correspond to a straight line distance to the next * volume intersection. In typical situations this leads to about - * one full leap-frog step to the next volume boundary. + * (at least) one full leap-frog step to the next volume boundary. * **/ -- GitLab