Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Calculate safety distances using OBZ when simple safety is not available #1427

Closed
wants to merge 2 commits into from
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 11 additions & 3 deletions src/orange/univ/SimpleUnitTracker.hh
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@
#include "corecel/math/Algorithms.hh"
#include "orange/OrangeData.hh"
#include "orange/detail/BIHTraverser.hh"
#include "orange/detail/OrientedBoundingZone.hh"
#include "orange/surf/LocalSurfaceVisitor.hh"

#include "detail/InfixEvaluator.hh"
Expand Down Expand Up @@ -100,6 +101,7 @@ class SimpleUnitTracker

ParamsRef const& params_;
SimpleUnitRecord const& unit_record_;
detail::OrientedBoundingZone::StoragePointers obz_storage_;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should probably be created on the fly since we keep params_ around


//// METHODS ////

Expand Down Expand Up @@ -143,7 +145,9 @@ class SimpleUnitTracker
*/
CELER_FUNCTION
SimpleUnitTracker::SimpleUnitTracker(ParamsRef const& params, SimpleUnitId suid)
: params_(params), unit_record_(params.simple_units[suid])
: params_(params)
, unit_record_(params.simple_units[suid])
, obz_storage_({&params_.transforms, &params_.reals})
{
CELER_EXPECT(params_);
}
Expand Down Expand Up @@ -301,8 +305,12 @@ CELER_FUNCTION real_type SimpleUnitTracker::safety(Real3 const& pos,
if (!vol.simple_safety())
{
// Has a tricky surface: we can't use the simple algorithm to calculate
// the safety, so return a conservative estimate.
return 0;
// the safety, use the oriented bounding zone
auto const& obz_id
= params_.volume_records[unit_record_.volumes[volid]].obz_id;
detail::OrientedBoundingZone obz(params_.obz_records[obz_id],
obz_storage_);
return obz.calc_safety_inside(pos);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Unfortunately it's not going to be this simple: we need the minimum over both "inside current volume" and "outside all directly enclosed volumes". I'm not sure if we have that metadata yet, in fact...

}

// Calculate minimim distance to all local faces
Expand Down
7 changes: 7 additions & 0 deletions test/TestMacros.hh
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,13 @@
# define TEST_IF_CELERITAS_GEANT(name) DISABLED_##name
#endif

//! Construct a test name that is disabled GEANT4 is disable or realtype=float
#if CELERITAS_USE_GEANT4 && CELERITAS_REAL_TYPE == CELERITAS_REAL_TYPE_DOUBLE
# define TEST_IF_CELERITAS_GEANT_AND_DOUBLE(name) name
#else
# define TEST_IF_CELERITAS_GEANT_AND_DOUBLE(name) DISABLED_##name
#endif

//! Construct a test name that is disabled when ROOT is disabled
#if CELERITAS_USE_ROOT
# define TEST_IF_CELERITAS_USE_ROOT(name) name
Expand Down
89 changes: 89 additions & 0 deletions test/geocel/data/cone.gdml
Original file line number Diff line number Diff line change
@@ -0,0 +1,89 @@
<?xml version="1.0"?>
<!-- Copyright 2020-2024 UT-Battelle, LLC, and other Celeritas developers. -->
<!-- See the top-level COPYRIGHT file for details. -->
<!-- SPDX-License-Identifier: (Apache-2.0 OR MIT) -->
<!-- \file cylinders.gdml -->
<!-- \brief three concentric cylinders -->
<gdml xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://service-spi.web.cern.ch/service-spi/app/releases/GDML/schema/gdml.xsd">

<materials>
<isotope N="1" Z="1" name="H1">
<atom unit="g/mole" value="1"/>
</isotope>
<element name="H">
<fraction n="1" ref="H1"/>
</element>
<material name="lo density celerogen" state="gas">
<T unit="K" value="2.73"/>
<P unit="pascal" value="3e-18"/>
<MEE unit="eV" value="21.8"/>
<D unit="g/cm3" value="1e-20"/>
<fraction n="1" ref="H"/>
</material>
<material name="hi density celerogen" state="gas">
<T unit="K" value="2.73"/>
<P unit="pascal" value="3e-18"/>
<MEE unit="eV" value="21.8"/>
<D unit="g/cm3" value="1e-19"/>
<fraction n="1" ref="H"/>
</material>
<material name="celer composite" state="solid">
<T unit="K" value="293.15"/>
<MEE unit="eV" value="173"/>
<D unit="g/cm3" value="1"/>
<fraction n="0.1" ref="H"/>
<fraction n="0.3" ref="C"/>
<fraction n="0.6" ref="O"/>
</material>
</materials>

<define>
<position name="origin" unit="cm" x="0" y="0" z="0" />
</define>

<solids>
<box name="OuterBox" lunit="cm" x="200" y="200" z="200"/>
<cone name="cone1" rmin1="0" rmax1="13" rmin2="0" rmax2="0" z="50" deltaphi="360" aunit="deg" lunit="cm"/>
<tube name="cyl2" rmax="20" z="50" deltaphi="360" aunit="deg" lunit="cm"/>
<tube name="cyl3" rmax="30" z="50" deltaphi="360" aunit="deg" lunit="cm"/>
</solids>

<structure>
<volume name="vol1">
<solidref ref="cone1"/>
<materialref ref="lo density celerogen"/>
</volume>
<volume name="vol2">
<solidref ref="cyl2"/>
<materialref ref="hi density celerogen"/>
</volume>
<volume name="vol3">
<solidref ref="cyl3"/>
<materialref ref="celer composite"/>
</volume>

<volume name="World">
<solidref ref="OuterBox"/>
<materialref ref="lo density celerogen"/>
<physvol>
<volumeref ref="vol1"/>
<position ref="origin"/>
</physvol>
<physvol>
<volumeref ref="vol2"/>
<position ref="origin"/>
</physvol>
<physvol>
<volumeref ref="vol3"/>
<position ref="origin"/>
</physvol>
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think these three enclosed volumes are overlapping. It might be more natural for you to make a geometry with orange/orangeinp/UnitProto.test.cc using the InputBuilderTest harness.

</volume>
</structure>

<setup name="Default" version="1.0">
<world ref="World"/>
</setup>


</gdml>

22 changes: 22 additions & 0 deletions test/orange/univ/SimpleUnitTracker.test.cc
Original file line number Diff line number Diff line change
Expand Up @@ -129,6 +129,12 @@
void SetUp() override { this->build_geometry("five-volumes.org.json"); }
};

#define ConeTest TEST_IF_CELERITAS_GEANT_AND_DOUBLE(ConeTest)
class ConeTest : public SimpleUnitTrackerTest
{
void SetUp() override { this->build_gdml_geometry("cone.gdml"); }
};

//---------------------------------------------------------------------------//
// TEST FIXTURE IMPLEMENTATION
//---------------------------------------------------------------------------//
Expand Down Expand Up @@ -967,6 +973,22 @@
}
}

TEST_F(ConeTest, safety)

Check failure on line 976 in test/orange/univ/SimpleUnitTracker.test.cc

View workflow job for this annotation

GitHub Actions / JUnit Test Report

/home/runner/work/celeritas/celeritas/test/orange/univ/SimpleUnitTracker.test.cc.safety

/home/runner/work/celeritas/celeritas/test/orange/univ/SimpleUnitTracker.test.cc:986 The difference between 13 and tracker.safety({0., 0., 0.}, vol1) is 117.00001525878906, which exceeds 1E-3, where 13 evaluates to 13, tracker.safety({0., 0., 0.}, vol1) evaluates to 130.00001525878906, and 1E-3 evaluates to 0.001.
Raw output
/home/runner/work/celeritas/celeritas/test/orange/univ/SimpleUnitTracker.test.cc:986
The difference between 13 and tracker.safety({0., 0., 0.}, vol1) is 117.00001525878906, which exceeds 1E-3, where
13 evaluates to 13,
tracker.safety({0., 0., 0.}, vol1) evaluates to 130.00001525878906, and
1E-3 evaluates to 0.001.
{
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh crap, I forgot that the Geant input gets scaled into the native unit system (unlike the .org.json files which are precalculated). This test is with the CLHEP unit system which scales everything by a factor of 10. I think it would be better to build a geometry with the InputBuilder test... sorry 😞

// Check safety distance within a cone, where simple safety is not
// supported
SimpleUnitTracker tracker(this->host_params(), SimpleUnitId{0});
detail::UniverseIndexer ui(this->host_params().universe_indexer_data);
LocalVolumeId vol1 = ui.local_volume(this->find_volume("vol1")).volume;

// In the center of the cone, the safety distance should be the radius of
// the cone
EXPECT_NEAR(13, tracker.safety({0., 0., 0.}, vol1), 1E-3);

// Outside the inner part of the OBZ, the safety distance should be zero
EXPECT_EQ(0., tracker.safety({0., 0., 1.}, vol1));
}

//---------------------------------------------------------------------------//
} // namespace test
} // namespace celeritas
Loading