IAP GITLAB

Commit 159919b1 authored by Ralf Ulrich's avatar Ralf Ulrich

fixed testCascade

parent d7c7484d
Pipeline #2710 failed with stages
in 15 minutes and 16 seconds
......@@ -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:
......
......@@ -19,6 +19,7 @@ set (
set (
ENVIRONMENT_HEADERS
VolumeTreeNode.h
IEmpty.hpp
IMediumModel.h
NuclearComposition.h
HomogeneousMedium.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>
......
......@@ -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;
......
......@@ -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;
......
......@@ -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);
}
......
......@@ -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>
......
......@@ -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;
}
......
......@@ -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
......
......@@ -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(
......
......@@ -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.
*
**/
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment