From 8ae37c8883bf90fba85301963cd27ba45820e5ba Mon Sep 17 00:00:00 2001
From: Christof Kraus <christof.kraus@psi.ch>
Date: Wed, 22 May 2019 19:48:58 +0200
Subject: [PATCH] Use different chains of random numbers on each core of a
parallel run in CollimatorPhysics. Furhtermore fix simulation time after all
particles were in material and fix spurious structure in longitudinal
position with period length of time step that was present after passage of
material.
---
src/Algorithms/ParallelTTracker.cpp | 32 +-
src/Algorithms/ParallelTTracker.h | 2 +-
src/Classic/AbsBeamline/Degrader.cpp | 7 +-
src/Classic/Algorithms/PartBunchBase.h | 3 +-
src/Classic/Algorithms/PartBunchBase.hpp | 6 +
src/Classic/Solvers/CollimatorPhysics.cpp | 795 ++++++++++--------
src/Classic/Solvers/CollimatorPhysics.hh | 42 +-
.../ParticleMatterInteractionHandler.hh | 4 +-
src/Structure/ParticleMatterInteraction.cpp | 23 +-
src/Structure/ParticleMatterInteraction.h | 13 +-
tools/emacs/opal.el | 4 +-
11 files changed, 495 insertions(+), 436 deletions(-)
diff --git a/src/Algorithms/ParallelTTracker.cpp b/src/Algorithms/ParallelTTracker.cpp
index 4d765acf3..f81e43be5 100644
--- a/src/Algorithms/ParallelTTracker.cpp
+++ b/src/Algorithms/ParallelTTracker.cpp
@@ -78,7 +78,7 @@ ParallelTTracker::ParallelTTracker(const Beamline &beamline,
fieldEvaluationTimer_m(IpplTimings::getTimer("External field eval")),
BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")),
WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
- particleMaterStatus_m(false),
+ particleMatterStatus_m(false),
totalParticlesInSimulation_m(0)
{
@@ -117,7 +117,7 @@ ParallelTTracker::ParallelTTracker(const Beamline &beamline,
fieldEvaluationTimer_m(IpplTimings::getTimer("External field eval")),
BinRepartTimer_m(IpplTimings::getTimer("Binaryrepart")),
WakeFieldTimer_m(IpplTimings::getTimer("WakeField")),
- particleMaterStatus_m(false),
+ particleMatterStatus_m(false),
totalParticlesInSimulation_m(0)
{
for (std::vector<unsigned long long>::const_iterator it = maxSteps.begin(); it != maxSteps.end(); ++ it) {
@@ -283,14 +283,14 @@ void ParallelTTracker::execute() {
setTime();
- double t = itsBunch_m->getT() - globalTimeShift;
- itsBunch_m->setT(t);
+ double time = itsBunch_m->getT() - globalTimeShift;
+ itsBunch_m->setT(time);
*gmsg << level1 << *itsBunch_m << endl;
unsigned long long step = itsBunch_m->getGlobalTrackStep();
OPALTimer::Timer myt1;
- *gmsg << "Track start at: " << myt1.time() << ", t= " << Util::getTimeString(t) << "; "
+ *gmsg << "Track start at: " << myt1.time() << ", t= " << Util::getTimeString(time) << "; "
<< "zstart at: " << Util::getLengthString(pathLength_m)
<< endl;
@@ -334,10 +334,9 @@ void ParallelTTracker::execute() {
timeIntegration2(pusher);
- t += itsBunch_m->getdT();
- itsBunch_m->setT(t);
+ itsBunch_m->incrementT();
- if (t > 0.0) {
+ if (itsBunch_m->getT() > 0.0) {
updateReference(pusher);
}
@@ -775,13 +774,13 @@ void ParallelTTracker::computeParticleMatterInteraction(IndexMap::value_t elemen
activeParticleMatterInteractionHandlers_m.insert(it);
}
- if(!particleMaterStatus_m) {
- msg << level2 << "============== START PARTICLE MATER INTERACTION CALCULATION =============" << endl;
- particleMaterStatus_m = true;
+ if(!particleMatterStatus_m) {
+ msg << level2 << "============== START PARTICLE MATTER INTERACTION CALCULATION =============" << endl;
+ particleMatterStatus_m = true;
}
}
- if (particleMaterStatus_m) {
+ if (particleMatterStatus_m) {
do {
///all particles in material if max per node is 2 and other degraders have 0 particles
//check if more than one degrader has particles
@@ -837,7 +836,8 @@ void ParallelTTracker::computeParticleMatterInteraction(IndexMap::value_t elemen
rediffusedParticles += it->getRediffused();
numEnteredParticles += it->getNumEntered();
- //if all particles where in material update time to time in degrader
+
+ //if all particles were in material update time to time in degrader
if (it->getFlagAllParticlesIn()) {
double timeDifference = it->getTime() - itsBunch_m->getdT() - itsBunch_m->getT();
if (timeDifference > 0.0) {
@@ -867,14 +867,14 @@ void ParallelTTracker::computeParticleMatterInteraction(IndexMap::value_t elemen
//if bunch has no particles update time to time in degrader
if (itsBunch_m->getTotalNum() == 0)
- itsBunch_m->setT(itsBunch_m->getT() + itsBunch_m->getdT());
+ itsBunch_m->incrementT();
} while (itsBunch_m->getTotalNum() == 0);
if (activeParticleMatterInteractionHandlers_m.size() == 0) {
- msg << level2 << "============== END PARTICLE MATER INTERACTION CALCULATION =============" << endl;
- particleMaterStatus_m = false;
+ msg << level2 << "============== END PARTICLE MATTER INTERACTION CALCULATION =============" << endl;
+ particleMatterStatus_m = false;
}
}
}
diff --git a/src/Algorithms/ParallelTTracker.h b/src/Algorithms/ParallelTTracker.h
index 3f8707bd7..dc35e2c71 100644
--- a/src/Algorithms/ParallelTTracker.h
+++ b/src/Algorithms/ParallelTTracker.h
@@ -246,7 +246,7 @@ private:
IpplTimings::TimerRef WakeFieldTimer_m;
std::set<ParticleMatterInteractionHandler*> activeParticleMatterInteractionHandlers_m;
- bool particleMaterStatus_m;
+ bool particleMatterStatus_m;
unsigned long totalParticlesInSimulation_m;
diff --git a/src/Classic/AbsBeamline/Degrader.cpp b/src/Classic/AbsBeamline/Degrader.cpp
index 2a6706552..7982966f8 100644
--- a/src/Classic/AbsBeamline/Degrader.cpp
+++ b/src/Classic/AbsBeamline/Degrader.cpp
@@ -131,10 +131,9 @@ bool Degrader::applyToReferenceParticle(const Vector_t &R,
Vector_t &B) {
if (!isInMaterial(R(2))) return false;
- const double eV2GeV = 1e-9;
- double Ekin = Util::getEnergy(P, RefPartBunch_m->getM() * eV2GeV);
- bool isDead = getParticleMatterInteraction()->computeEnergyLoss(Ekin, RefPartBunch_m->getdT(), false);
- double deltaP = Util::getP(Ekin, RefPartBunch_m->getM() * eV2GeV) - euclidean_norm(P);
+ Vector_t updatedP = P;
+ bool isDead = getParticleMatterInteraction()->computeEnergyLoss(updatedP, RefPartBunch_m->getdT(), false);
+ double deltaP = euclidean_norm(updatedP) - euclidean_norm(P);
E(2) += deltaP * RefPartBunch_m->getM() / (RefPartBunch_m->getdT() * RefPartBunch_m->getQ() * Physics::c);
return isDead;
diff --git a/src/Classic/Algorithms/PartBunchBase.h b/src/Classic/Algorithms/PartBunchBase.h
index d7212ee38..fd6933810 100644
--- a/src/Classic/Algorithms/PartBunchBase.h
+++ b/src/Classic/Algorithms/PartBunchBase.h
@@ -208,6 +208,7 @@ public:
double getdT() const;
void setT(double t);
+ void incrementT();
double getT() const;
/**
@@ -675,4 +676,4 @@ protected:
#include "PartBunchBase.hpp"
-#endif
+#endif
\ No newline at end of file
diff --git a/src/Classic/Algorithms/PartBunchBase.hpp b/src/Classic/Algorithms/PartBunchBase.hpp
index d8ae6f9a4..4e296cc02 100644
--- a/src/Classic/Algorithms/PartBunchBase.hpp
+++ b/src/Classic/Algorithms/PartBunchBase.hpp
@@ -1076,6 +1076,12 @@ void PartBunchBase<T, Dim>::setT(double t) {
}
+template <class T, unsigned Dim>
+void PartBunchBase<T, Dim>::incrementT() {
+ t_m += dt_m;
+}
+
+
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getT() const {
return t_m;
diff --git a/src/Classic/Solvers/CollimatorPhysics.cpp b/src/Classic/Solvers/CollimatorPhysics.cpp
index d41227396..3601fe62a 100644
--- a/src/Classic/Solvers/CollimatorPhysics.cpp
+++ b/src/Classic/Solvers/CollimatorPhysics.cpp
@@ -17,7 +17,7 @@
#include "AbsBeamline/SBend.h"
#include "AbsBeamline/RBend.h"
#include "AbsBeamline/Multipole.h"
-#include "Structure/LossDataSink.h" // OPAL file
+#include "Structure/LossDataSink.h"
#include "Utilities/Options.h"
#include "Utilities/GeneralClassicException.h"
#include "Utilities/Util.h"
@@ -30,13 +30,6 @@
#include <algorithm>
using Physics::pi;
-using Physics::m_p;
-using Physics::m_e;
-using Physics::h_bar;
-using Physics::epsilon_0;
-using Physics::r_e;
-using Physics::z_p;
-using Physics::Avo;
#ifdef OPAL_DKS
const int CollimatorPhysics::numpar = 13;
@@ -83,7 +76,10 @@ namespace {
};
}
-CollimatorPhysics::CollimatorPhysics(const std::string &name, ElementBase *element, std::string &material):
+CollimatorPhysics::CollimatorPhysics(const std::string &name,
+ ElementBase *element,
+ std::string &material,
+ bool enableRutherford):
ParticleMatterInteractionHandler(name, element),
T_m(0.0),
dT_m(0.0),
@@ -101,10 +97,11 @@ CollimatorPhysics::CollimatorPhysics(const std::string &name, ElementBase *eleme
bunchToMatStat_m(0),
stoppedPartStat_m(0),
rediffusedStat_m(0),
- locPartsInMat_m(0),
+ totalPartsInMat_m(0),
Eavg_m(0.0),
Emax_m(0.0),
- Emin_m(0.0)
+ Emin_m(0.0),
+ enableRutherford_m(enableRutherford)
#ifdef OPAL_DKS
, curandInitSet(0)
, ierr(0)
@@ -118,7 +115,16 @@ CollimatorPhysics::CollimatorPhysics(const std::string &name, ElementBase *eleme
gsl_rng_env_setup();
rGen_m = gsl_rng_alloc(gsl_rng_default);
- gsl_rng_set(rGen_m, Options::seed);
+
+ size_t mySeed = Options::seed;
+
+ if (Options::seed == -1) {
+ struct timeval tv;
+ gettimeofday(&tv,0);
+ mySeed = tv.tv_sec + tv.tv_usec;
+ }
+
+ gsl_rng_set(rGen_m, mySeed + Ippl::myNode());
configureMaterialParameters();
@@ -172,115 +178,22 @@ std::string CollimatorPhysics::getName() {
return element_ref_m->getName();
}
-void CollimatorPhysics::doPhysics(PartBunchBase<double, 3> *bunch) {
- /***
- Do physics if
- -- particle in material
- -- particle not dead (locParts_m[i].label != -1.0)
-
- Absorbed particle i: locParts_m[i].label = -1.0;
-
- Particle goes back to beam if
- -- not absorbed and out of material
- */
-
- for (size_t i = 0; i < locParts_m.size(); ++i) {
- Vector_t &R = locParts_m[i].Rincol;
- Vector_t &P = locParts_m[i].Pincol;
-
- double Eng = (sqrt(1.0 + dot(P, P)) - 1) * m_p;
- if (locParts_m[i].label != -1) {
-
- if (hitTester_m->checkHit(R, P, dT_m)) {
- bool pdead = computeEnergyLoss(Eng, dT_m);
- if (!pdead) {
- double ptot = sqrt((m_p + Eng) * (m_p + Eng) - (m_p) * (m_p)) / m_p;
- P = ptot * P / sqrt(dot(P, P));
- /*
- Now scatter and transport particle in material.
- The checkHit call just above will detect if the
- particle is rediffused from the material into vacuum.
- */
-
- computeCoulombScattering(R, P, dT_m);
- } else {
- // The particle is stopped in the material, set label to -1
- locParts_m[i].label = -1.0;
- stoppedPartStat_m++;
- lossDs_m->addParticle(R, P, -locParts_m[i].IDincol);
- }
- } else {
- /* The particle exits the material but is still in the loop of the substep,
- Finish the timestep by letting the particle drift and after the last
- substep call addBackToBunch
- */
- double gamma = (Eng + m_p) / m_p;
- double beta = sqrt(1.0 - 1.0 / (gamma * gamma));
- if (collshape_m == ElementBase::CCOLLIMATOR) {
- R = R + dT_m * beta * Physics::c * P / sqrt(dot(P, P)) * 1000;
- } else {
- R = R + dT_m * Physics::c * P / sqrt(1.0 + dot(P, P)) ;
- addBackToBunch(bunch, i);
- }
- }
- }
- }
-}
-
-/// Energy Loss: using the Bethe-Bloch equation.
-/// Energy straggling: For relatively thick absorbers such that the number of collisions is large,
-/// the energy loss distribution is shown to be Gaussian in form.
-/// See Particle Physics Booklet, chapter 'Passage of particles through matter'.
-// -------------------------------------------------------------------------
-bool CollimatorPhysics::computeEnergyLoss(double &Eng, const double deltat, bool includeFluctuations) const {
- /// Eng GeV
-
- const double m2cm = 1e2;
- const double GeV2keV = 1e6;
-
- Eng *= GeV2keV;
- double dEdx = 0.0;
-
- const double gamma = (Eng + m_p * GeV2keV) / (m_p * GeV2keV);
- const double beta = sqrt(1.0 - 1.0 / (gamma * gamma));
-
- const double deltas = deltat * beta * Physics::c;
- const double deltasrho = deltas * m2cm * rho_m;
- const double K = 4.0 * pi * Avo * std::pow(r_e * m2cm, 2) * m_e * GeV2keV;
- const double sigma_E = sqrt(K * m_e * GeV2keV * rho_m * (Z_m / A_m) * deltas * m2cm);
-
- if (Eng >= 600) {
- const double eV2keV = 1e-3;
- const double gammaSqr = gamma * gamma;
- const double betaSqr = beta * beta;
- const double Tmax = (2.0 * m_e * GeV2keV * betaSqr * gammaSqr /
- (std::pow(gamma + m_e / m_p, 2) - (gammaSqr - 1.0)));
-
- dEdx = (-K * z_p * z_p * Z_m / (A_m * betaSqr) *
- (1.0 / 2.0 * std::log(2 * m_e * GeV2keV * betaSqr * gammaSqr * Tmax / std::pow(I_m * eV2keV, 2)) - betaSqr));
-
- Eng += deltasrho * dEdx;
- } else if (Eng > 10) {
- const double Ts = Eng / 1.0073; // 1.0073 is the mass of proton in atomic mass units.
- const double epsilon_low = A2_c * pow(Ts, 0.45);
- const double epsilon_high = (A3_c / Ts) * log(1 + (A4_c / Ts) + (A5_c * Ts));
- const double epsilon = (epsilon_low * epsilon_high) / (epsilon_low + epsilon_high);
- dEdx = -epsilon / (1e18 * (A_m / Avo));
-
- Eng += deltasrho * dEdx;
- }
-
- if (includeFluctuations)
- Eng += gsl_ran_gaussian(rGen_m, sigma_E);
-
- bool dead = (Eng < 10 || dEdx > 0);
-
- Eng /= GeV2keV;
+/// The material of the collimator
+// ------------------------------------------------------------------------
+void CollimatorPhysics::configureMaterialParameters() {
- return dead;
+ auto material = Physics::Material::getMaterial(material_m);
+ Z_m = material->getAtomicNumber();
+ A_m = material->getAtomicMass();
+ rho_m = material->getMassDensity();
+ X0_m = material->getRadiationLength();
+ I_m = material->getMeanExcitationEnergy();
+ A2_c = material->getStoppingPowerFitCoefficients(Physics::Material::A2);
+ A3_c = material->getStoppingPowerFitCoefficients(Physics::Material::A3);
+ A4_c = material->getStoppingPowerFitCoefficients(Physics::Material::A4);
+ A5_c = material->getStoppingPowerFitCoefficients(Physics::Material::A5);
}
-
void CollimatorPhysics::apply(PartBunchBase<double, 3> *bunch,
const std::pair<Vector_t, double> &boundingSphere,
size_t numParticlesInSimulation) {
@@ -288,15 +201,11 @@ void CollimatorPhysics::apply(PartBunchBase<double, 3> *bunch,
/*
Particles that have entered material are flagged as Bin[i] == -1.
- Fixme: should use PType
Flagged particles are copied to a local structure within Collimator Physics locParts_m.
Particles in that structure will be pushed in the material and either come
- back to the bunch or will be fully stopped in the material. For the push in the
- material we use sub-timesteps.
-
- Newly entered particles will be copied to locParts_m at the end of apply.
+ back to the bunch or will be fully stopped in the material.
*/
Eavg_m = 0.0;
@@ -306,15 +215,11 @@ void CollimatorPhysics::apply(PartBunchBase<double, 3> *bunch,
bunchToMatStat_m = 0;
rediffusedStat_m = 0;
stoppedPartStat_m = 0;
- locPartsInMat_m = 0;
+ totalPartsInMat_m = 0;
dT_m = bunch->getdT();
T_m = bunch->getT();
- /*
- Because this is not properly set in the Component class when calling in the Constructor
- */
-
#ifdef OPAL_DKS
if (collshape_m == ElementBase::DEGRADER && IpplInfo::DKSEnabled) {
applyDKS(bunch, boundingSphere, numParticlesInSimulation);
@@ -328,282 +233,282 @@ void CollimatorPhysics::apply(PartBunchBase<double, 3> *bunch,
IpplTimings::stopTimer(DegraderApplyTimer_m);
}
-#ifdef OPAL_DKS
-void CollimatorPhysics::applyDKS(PartBunchBase<double, 3> *bunch,
- const std::pair<Vector_t, double> &boundingSphere,
- size_t numParticlesInSimulation) {
+void CollimatorPhysics::applyNonDKS(PartBunchBase<double, 3> *bunch,
+ const std::pair<Vector_t, double> &boundingSphere,
+ size_t numParticlesInSimulation) {
bool onlyOneLoopOverParticles = ! (allParticleInMat_m);
- //if firs call to apply setup needed accelerator resources
- setupCollimatorDKS(bunch, numParticlesInSimulation);
-
do {
IpplTimings::startTimer(DegraderLoopTimer_m);
+ push();
+ setTimeStepForLeavingParticles();
- //write particles to GPU if there are any to write
- if (dksParts_m.size() > 0) {
- //wrtie data from dksParts_m to the end of mem_ptr (offset = numparticles)
- dksbase.writeDataAsync<PART_DKS>(mem_ptr, &dksParts_m[0],
- dksParts_m.size(), -1, numparticles);
+ /*
+ if we are not looping copy newly arrived particles
+ */
+ if (onlyOneLoopOverParticles)
+ copyFromBunch(bunch, boundingSphere);
- //update number of particles on Device
- numparticles += dksParts_m.size();
+ addBackToBunch(bunch);
- //free locParts_m vector
- dksParts_m.erase(dksParts_m.begin(), dksParts_m.end());
- }
+ computeInteraction();
- //execute CollimatorPhysics kernels on GPU if any particles are there
- if (numparticles > 0) {
- dksbase.callCollimatorPhysics2(mem_ptr, par_ptr, numparticles);
- }
+ push();
+ resetTimeStep();
- //sort device particles and get number of particles comming back to bunch
- int numaddback = 0;
- if (numparticles > 0) {
- dksbase.callCollimatorPhysicsSort(mem_ptr, numparticles, numaddback);
+ IpplTimings::stopTimer(DegraderLoopTimer_m);
+
+ T_m += dT_m; // update local time
+
+ unsigned int locPartsInMat = locParts_m.size();
+ reduce(locPartsInMat, totalPartsInMat_m, OpAddAssign());
+
+ if (allParticleInMat_m) {
+ onlyOneLoopOverParticles = (rediffusedStat_m > 0 || totalPartsInMat_m <= 0);
+ } else {
+ onlyOneLoopOverParticles = true;
}
+ } while (onlyOneLoopOverParticles == false);
+}
- //read particles from GPU if any are comming out of material
- if (numaddback > 0) {
+void CollimatorPhysics::computeInteraction() {
+ /***
+ Do physics if
+ -- particle not stopped (locParts_m[i].label != -1.0)
- //resize dksParts_m to hold particles that need to go back to bunch
- dksParts_m.resize(numaddback);
+ Absorbed particle i: locParts_m[i].label = -1.0;
+ */
- //read particles that need to be added back to bunch
- //particles that need to be added back are at the end of Device array
- dksbase.readData<PART_DKS>(mem_ptr, &dksParts_m[0], numaddback,
- numparticles - numaddback);
+ for (size_t i = 0; i < locParts_m.size(); ++i) {
+ if (locParts_m[i].label != -1) {
+ Vector_t &R = locParts_m[i].Rincol;
+ Vector_t &P = locParts_m[i].Pincol;
+ double &dt = locParts_m[i].DTincol;
- //add particles back to the bunch
- for (unsigned int i = 0; i < dksParts_m.size(); ++i) {
- if (dksParts_m[i].label == -2) {
- addBackToBunchDKS(bunch, i);
- rediffusedStat_m++;
+ if (hitTester_m->checkHit(R, P, dt)) {
+ bool pdead = computeEnergyLoss(P, dt);
+ if (!pdead) {
+ /*
+ Now scatter and transport particle in material.
+ The checkHit call just above will detect if the
+ particle is rediffused from the material into vacuum.
+ */
+
+ computeCoulombScattering(R, P, dt);
} else {
- stoppedPartStat_m++;
- lossDs_m->addParticle(dksParts_m[i].Rincol, dksParts_m[i].Pincol,
- -locParts_m[dksParts_m[i].localID].IDincol);
+ // The particle is stopped in the material, set label to -1
+ locParts_m[i].label = -1.0;
+ ++ stoppedPartStat_m;
+ lossDs_m->addParticle(R, P, -locParts_m[i].IDincol);
}
}
-
- //erase particles that came from device from host array
- dksParts_m.erase(dksParts_m.begin(), dksParts_m.end());
-
- //update number of particles on Device
- numparticles -= numaddback;
}
+ }
- IpplTimings::stopTimer(DegraderLoopTimer_m);
-
- if (onlyOneLoopOverParticles)
- copyFromBunchDKS(bunch, boundingSphere);
+ /*
+ delete absorbed particles
+ */
+ deleteParticleFromLocalVector();
+}
- T_m += dT_m;
+/// Energy Loss: using the Bethe-Bloch equation.
+/// Energy straggling: For relatively thick absorbers such that the number of collisions is large,
+/// the energy loss distribution is shown to be Gaussian in form.
+/// See Particle Physics Booklet, chapter 'Passage of particles through matter' or
+/// Review of Particle Physics, DOI: 10.1103/PhysRevD.86.010001, page 329 ff
+// -------------------------------------------------------------------------
+bool CollimatorPhysics::computeEnergyLoss(Vector_t &P,
+ const double deltat,
+ bool includeFluctuations) const {
+
+ constexpr double m2cm = 1e2;
+ constexpr double GeV2keV = 1e6;
+ constexpr double massElectron_keV = Physics::m_e * GeV2keV;
+ constexpr double massProton_keV = Physics::m_p * GeV2keV;
+ constexpr double massProton_amu = 1.007276466879;
+ constexpr double chargeProton = Physics::z_p;
+ constexpr double K = 4.0 * pi * Physics::Avo * std::pow(Physics::r_e * m2cm, 2) * massElectron_keV;
+
+ double gamma = Util::getGamma(P);
+ const double gammaSqr = std::pow(gamma, 2);
+ const double betaSqr = 1.0 - 1.0 / gammaSqr;
+ double beta = sqrt(betaSqr);
+ double Ekin = (gamma - 1) * massProton_keV;
+ double dEdx = 0.0;
- locPartsInMat_m = numparticles;
- reduce(locPartsInMat_m, locPartsInMat_m, OpAddAssign());
+ const double deltas = deltat * beta * Physics::c;
+ const double deltasrho = deltas * m2cm * rho_m;
- int maxPerNode = bunch->getLocalNum();
- reduce(maxPerNode, maxPerNode, OpMaxAssign());
+ if (Ekin >= 600) {
+ constexpr double eV2keV = 1e-3;
+ constexpr double massRatio = Physics::m_e / Physics::m_p;
+ double Tmax = (2.0 * massElectron_keV * betaSqr * gammaSqr /
+ (std::pow(gamma + massRatio, 2) - (gammaSqr - 1.0)));
- //more than one loop only if all the particles are in this degrader
- if (allParticleInMat_m) {
- onlyOneLoopOverParticles = ( (unsigned)maxPerNode > bunch->getMinimumNumberOfParticlesPerCore() || locPartsInMat_m <= 0);
- } else {
- onlyOneLoopOverParticles = true;
- }
+ dEdx = (-K * std::pow(chargeProton, 2) * Z_m / (A_m * betaSqr) *
+ (0.5 * std::log(2 * massElectron_keV * betaSqr * gammaSqr * Tmax / std::pow(I_m * eV2keV, 2)) - betaSqr));
- } while (onlyOneLoopOverParticles == false);
-}
-#endif
+ Ekin += deltasrho * dEdx;
+ } else if (Ekin > 10) {
+ const double Ts = Ekin / massProton_amu;
+ const double epsilon_low = A2_c * pow(Ts, 0.45);
+ const double epsilon_high = (A3_c / Ts) * log(1 + (A4_c / Ts) + (A5_c * Ts));
+ const double epsilon = (epsilon_low * epsilon_high) / (epsilon_low + epsilon_high);
+ dEdx = -epsilon / (1e18 * (A_m / Physics::Avo));
-void CollimatorPhysics::applyNonDKS(PartBunchBase<double, 3> *bunch,
- const std::pair<Vector_t, double> &boundingSphere,
- size_t numParticlesInSimulation) {
- bool onlyOneLoopOverParticles = ! (allParticleInMat_m);
+ Ekin += deltasrho * dEdx;
+ }
- do {
- IpplTimings::startTimer(DegraderLoopTimer_m);
- doPhysics(bunch);
- /*
- delete absorbed particles and particles that went to the bunch
- */
- deleteParticleFromLocalVector();
- IpplTimings::stopTimer(DegraderLoopTimer_m);
+ if (includeFluctuations) {
+ double sigma_E = sqrt(K * massElectron_keV * rho_m * (Z_m / A_m) * deltas * m2cm);
+ Ekin += gsl_ran_gaussian(rGen_m, sigma_E);
+ }
- /*
- if we are not looping copy newly arrived particles
- */
- if (onlyOneLoopOverParticles)
- copyFromBunch(bunch, boundingSphere);
+ gamma = Ekin / massProton_keV + 1.0;
+ beta = sqrt(1.0 - 1.0 / std::pow(gamma, 2));
+ P = gamma * beta * P / euclidean_norm(P);
- T_m += dT_m; // update local time
+ bool stopped = (Ekin < 10 || dEdx > 0);
+ return stopped;
+}
- locPartsInMat_m = locParts_m.size();
- reduce(locPartsInMat_m, locPartsInMat_m, OpAddAssign());
- int maxPerNode = bunch->getLocalNum();
- reduce(maxPerNode, maxPerNode, OpMaxAssign());
- if (allParticleInMat_m) {
- onlyOneLoopOverParticles = ( (unsigned)maxPerNode > bunch->getMinimumNumberOfParticlesPerCore() || locPartsInMat_m <= 0);
- } else {
- onlyOneLoopOverParticles = true;
- }
- } while (onlyOneLoopOverParticles == false);
+// Implement the rotation in 2 dimensions here
+// For details see:
+// J. Beringer et al. (Particle Data Group), "Passage of particles through matter"
+// Phys. Rev. D 86, 010001 (2012),
+void CollimatorPhysics::applyRotation(Vector_t &P,
+ Vector_t &R,
+ double shift,
+ double thetacou) {
+ // Calculate the angle between the transverse and longitudinal component of the momentum
+ double Psixz = std::fmod(std::atan2(P(0), P(2)) + Physics::two_pi, Physics::two_pi);
+
+ R(0) = R(0) + shift * cos(Psixz);
+ R(2) = R(2) - shift * sin(Psixz);
+
+ // Apply the rotation about the random angle thetacou
+ double Px = P(0);
+ P(0) = Px * cos(thetacou) + P(2) * sin(thetacou);
+ P(2) = -Px * sin(thetacou) + P(2) * cos(thetacou);
}
-/// The material of the collimator
-// ------------------------------------------------------------------------
-void CollimatorPhysics::configureMaterialParameters() {
+void CollimatorPhysics::applyRandomRotation(Vector_t &P, double theta0) {
- auto material = Physics::Material::getMaterial(material_m);
- Z_m = material->getAtomicNumber();
- A_m = material->getAtomicMass();
- rho_m = material->getMassDensity();
- X0_m = material->getRadiationLength();
- I_m = material->getMeanExcitationEnergy();
- A2_c = material->getStoppingPowerFitCoefficients(Physics::Material::A2);
- A3_c = material->getStoppingPowerFitCoefficients(Physics::Material::A3);
- A4_c = material->getStoppingPowerFitCoefficients(Physics::Material::A4);
- A5_c = material->getStoppingPowerFitCoefficients(Physics::Material::A5);
-}
+ double thetaru = 2.5 / sqrt(gsl_rng_uniform(rGen_m)) * 2.0 * theta0;
+ double phiru = Physics::two_pi * gsl_rng_uniform(rGen_m);
-// Implement the rotation in 2 dimensions here
-// For details: see J. Beringer et al. (Particle Data Group), Phys. Rev. D 86, 010001 (2012), "Passage of particles through matter"
-void CollimatorPhysics::applyRotation(double &px,
- double &pz,
- double &x,
- double &z,
- double xplane,
- double normP,
- double thetacou,
- double deltas,
- int coord) {
- // Calculate the angle between the px and pz momenta to change from beam coordinate to lab coordinate
- double Psixz = std::fmod(std::atan2(px,pz) + Physics::two_pi, Physics::two_pi);
-
- // Apply the rotation about the random angle thetacou & change from beam
- // coordinate system to the lab coordinate system using Psixz (2 dimensions)
- if (coord == 1) {
- x = x + deltas * px / normP + xplane * cos(Psixz);
- z = z - xplane * sin(Psixz);
- }
- if (coord == 2) {
- x = x + deltas * px / normP + xplane * cos(Psixz);
- z = z - xplane * sin(Psixz) + deltas * pz / normP;
- }
+ double normPtrans = sqrt(P(0) * P(0) + P(1) * P(1));
+ double Theta = std::atan(normPtrans / std::abs(P(2)));
+ double CosT = cos(Theta);
+ double SinT = sin(Theta);
- double pxz = sqrt(px * px + pz * pz);
- px = pxz * cos(Psixz) * sin(thetacou) + pxz * sin(Psixz) * cos(thetacou);
- pz = -pxz * sin(Psixz) * sin(thetacou) + pxz * cos(Psixz) * cos(thetacou);
-}
+ Vector_t X(cos(phiru)*sin(thetaru),
+ sin(phiru)*sin(thetaru),
+ cos(thetaru));
+ X *= euclidean_norm(P);
+
+ Vector_t W(-P(1), P(0), 0.0);
+ W = W / normPtrans;
-Vector_t ArbitraryRotation(Vector_t &W, Vector_t &Rorg, double Theta) {
- double C=cos(Theta);
- double S=sin(Theta);
- W = W / sqrt(dot(W,W));
- return Rorg * C + cross(W,Rorg) * S + W * dot(W,Rorg) * (1.0-C);
+ // Rodrigues' formula for a rotation about W by angle Theta
+ P = X * CosT + cross(W, X) * SinT + W * dot(W, X) * (1.0 - CosT);
}
/// Coulomb Scattering: Including Multiple Coulomb Scattering and large angle Rutherford Scattering.
/// Using the distribution given in Classical Electrodynamics, by J. D. Jackson.
//--------------------------------------------------------------------------
-void CollimatorPhysics::computeCoulombScattering(Vector_t &R, Vector_t &P, const double &deltat) {
+void CollimatorPhysics::computeCoulombScattering(Vector_t &R,
+ Vector_t &P,
+ double dt) {
+
+ constexpr double GeV2eV = 1e9;
+ constexpr double massProton_eV = Physics::m_p * GeV2eV;
+ constexpr double chargeProton = Physics::z_p;
+ constexpr double sqrtThreeInv = 1.0 / sqrt(3.0);
+ const double normP = euclidean_norm(P);
+ const double gammaSqr = std::pow(normP, 2) + 1.0;
+ const double beta = sqrt(1.0 - 1.0 / gammaSqr);
+ const double deltas = dt * beta * Physics::c;
+ const double theta0 = (13.6e6 / (beta * normP * massProton_eV) *
+ chargeProton * sqrt(deltas / X0_m) *
+ (1.0 + 0.038 * log(deltas / X0_m)));
+
+ double phi = Physics::two_pi * gsl_rng_uniform(rGen_m);
+ for (unsigned int i = 0; i < 2; ++ i) {
+ CoordinateSystemTrafo randomTrafo(R, Quaternion(cos(phi), 0, 0, sin(phi)));
+ P = randomTrafo.rotateTo(P);
+ R = Vector_t(0.0); // corresponds to randomTrafo.transformTo(R);
+
+ double z1 = gsl_ran_gaussian(rGen_m, 1.0);
+ double z2 = gsl_ran_gaussian(rGen_m, 1.0);
+
+ while(std::abs(z2) > 3.5) {
+ z1 = gsl_ran_gaussian(rGen_m, 1.0);
+ z2 = gsl_ran_gaussian(rGen_m, 1.0);
+ }
- const double Eng = sqrt(dot(P, P) + 1.0) * m_p - m_p;
- const double gamma = (Eng + m_p) / m_p;
- const double beta = sqrt(1.0 - 1.0 / (gamma * gamma));
- const double normP = sqrt(dot(P, P));
- const double deltas = deltat * beta * Physics::c;
- const double theta0 = 13.6e6 / (beta * normP * m_p * 1e9) * z_p * sqrt(deltas / X0_m) * (1.0 + 0.038 * log(deltas / X0_m));
+ double thetacou = z2 * theta0;
+ double shift = (z1 * sqrtThreeInv + z2) * deltas * theta0 * 0.5;
+ applyRotation(P, R, shift, thetacou);
- // x-direction: See Physical Review, "Multiple Scattering"
- double z1 = gsl_ran_gaussian(rGen_m, 1.0);
- double z2 = gsl_ran_gaussian(rGen_m, 1.0);
- double thetacou = z2 * theta0;
+ P = randomTrafo.rotateFrom(P);
+ R = randomTrafo.transformFrom(R);
- while(fabs(thetacou) > 3.5 * theta0) {
- z1 = gsl_ran_gaussian(rGen_m, 1.0);
- z2 = gsl_ran_gaussian(rGen_m, 1.0);
- thetacou = z2 * theta0;
- }
- double xplane = (z1 / sqrt(3.0) + z2) * deltas * theta0 / 2.0;
- // Apply change in coordinates for multiple scattering but not for Rutherford
- // scattering (take the deltas step only once each turn)
- int coord = 1;
- applyRotation(P(0), P(2), R(0), R(2), xplane, normP, thetacou, deltas, coord);
-
-
- // Rutherford-scattering in x-direction
- if (collshape_m == ElementBase::CCOLLIMATOR)
- R = R * 1e-3; // mm -> m
-
- // y-direction: See Physical Review, "Multiple Scattering"
- z1 = gsl_ran_gaussian(rGen_m, 1.0);
- z2 = gsl_ran_gaussian(rGen_m, 1.0);
- thetacou = z2 * theta0;
-
- while(fabs(thetacou) > 3.5 * theta0) {
- z1 = gsl_ran_gaussian(rGen_m, 1.0);
- z2 = gsl_ran_gaussian(rGen_m, 1.0);
- thetacou = z2 * theta0;
+ phi += 0.5 * Physics::pi;
}
- double yplane = z1 * deltas * theta0 / sqrt(12.0) + z2 * deltas * theta0 / 2.0;
- // Apply change in coordinates for multiple scattering but not for Rutherford
- // scattering (take the deltas step only once each turn)
- coord = 2;
- applyRotation(P(1), P(2), R(1), R(2), yplane, normP, thetacou, deltas, coord);
-
- // Rutherford-scattering in x-direction
- if (collshape_m == ElementBase::CCOLLIMATOR)
- R = R * 1e3; // m -> mm
-
- double P2 = gsl_rng_uniform(rGen_m);
- if (P2 < 0.0047) {
- double P3 = gsl_rng_uniform(rGen_m);
-
- double thetaru = 2.5 * sqrt(1 / P3) * 2.0 * theta0;
- double phiru = 2.0 * M_PI * gsl_rng_uniform(rGen_m);
- double th0=atan2(sqrt(P(0)*P(0)+P(1)*P(1)),fabs(P(2)));
- Vector_t W,X;
- X(0)=cos(phiru)*sin(thetaru);
- X(1)=sin(phiru)*sin(thetaru);
- X(2)=cos(thetaru);
- X*=sqrt(dot(P,P));
- W(0)=-P(1);
- W(1)= P(0);
- W(2)= 0.0;
- P = ArbitraryRotation(W, X,th0);
+
+ if (enableRutherford_m &&
+ gsl_rng_uniform(rGen_m) < 0.0047) {
+
+ applyRandomRotation(P, theta0);
}
}
-void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3> *bunch, unsigned i) {
+void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3> *bunch) {
- bunch->createWithID(locParts_m[i].IDincol);
+ const size_t nL = locParts_m.size();
+ if (nL == 0) return;
- /*
- Binincol is still <0, but now the particle is rediffused
- from the material and hence this is not a "lost" particle anymore
- */
+ unsigned int numLocalParticles = bunch->getLocalNum();
+ const double elementLength = element_ref_m->getElementLength();
- const unsigned int last = bunch->getLocalNum() - 1;
- bunch->Bin[last] = 1;
- bunch->R[last] = locParts_m[i].Rincol;
- bunch->P[last] = locParts_m[i].Pincol;
- bunch->Q[last] = locParts_m[i].Qincol;
- bunch->Bf[last] = locParts_m[i].Bfincol;
- bunch->Ef[last] = locParts_m[i].Efincol;
- bunch->dt[last] = locParts_m[i].DTincol;
+ for (size_t i = 0; i < nL; ++ i) {
+ Vector_t &R = locParts_m[i].Rincol;
+
+ if (R[2] >= elementLength) {
+
+ bunch->createWithID(locParts_m[i].IDincol);
+
+ /*
+ Binincol is still <0, but now the particle is rediffused
+ from the material and hence this is not a "lost" particle anymore
+ */
+
+ bunch->Bin[numLocalParticles] = 1;
+ bunch->R[numLocalParticles] = R;
+ bunch->P[numLocalParticles] = locParts_m[i].Pincol;
+ bunch->Q[numLocalParticles] = locParts_m[i].Qincol;
+ bunch->Bf[numLocalParticles] = 0.0;
+ bunch->Ef[numLocalParticles] = 0.0;
+ bunch->dt[numLocalParticles] = dT_m;
+
+ /*
+ This particle is back to the bunch, by set
+ ting the label to -1.0
+ the particle will be deleted.
+ */
+ locParts_m[i].label = -1.0;
+
+ ++ rediffusedStat_m;
+ ++ numLocalParticles;
+ }
+ }
/*
- This particle is back to the bunch, by set
- ting the label to -1.0
- the particle will be deleted.
+ delete particles that went to the bunch
*/
- locParts_m[i].label = -1.0;
-
- ++ rediffusedStat_m;
+ deleteParticleFromLocalVector();
}
void CollimatorPhysics::copyFromBunch(PartBunchBase<double, 3> *bunch,
@@ -628,9 +533,23 @@ void CollimatorPhysics::copyFromBunch(PartBunchBase<double, 3> *bunch,
((nL - ne) > minNumOfParticlesPerCore) &&
hitTester_m->checkHit(bunch->R[i], bunch->P[i], dT_m))
{
+ // adjust the time step for those particles that enter the material
+ // such that it corresponds to the time needed to reach the curren
+ // location form the edge of the material. Only use this time step
+ // for the computation of the interaction with the material, not for
+ // the integration of the particles. This will ensure that the momenta
+ // of all particles are reduced by approcimately the same amount in
+ // computeEnergyLoss.
+ double betaz = bunch->P[i](2) / Util::getGamma(bunch->P[i]);
+ double stepWidth = betaz * Physics::c * bunch->dt[i];
+ double tau = bunch->R[i](2) / stepWidth;
+
+ PAssert_LE(tau, 1.0);
+ PAssert_GE(tau, 0.0);
+
PART x;
x.localID = i;
- x.DTincol = bunch->dt[i];
+ x.DTincol = bunch->dt[i] * tau;
x.IDincol = bunch->ID[i];
x.Binincol = bunch->Bin[i];
x.Rincol = bunch->R[i];
@@ -663,20 +582,33 @@ void CollimatorPhysics::print(Inform &msg) {
Inform::FmtFlags_t ff = msg.flags();
// ToDo: need to move that to a statistics function
+ unsigned int locPartsInMat;
#ifdef OPAL_DKS
if (collshape_m == ElementBase::DEGRADER && IpplInfo::DKSEnabled)
locPartsInMat_m = numparticles + dksParts_m.size();
else
- locPartsInMat_m = locParts_m.size();
+ locPartsInMat = locParts_m.size();
#else
- locPartsInMat_m = locParts_m.size();
+ locPartsInMat = locParts_m.size();
#endif
- reduce(locPartsInMat_m, locPartsInMat_m, OpAddAssign());
- reduce(bunchToMatStat_m, bunchToMatStat_m, OpAddAssign());
- reduce(rediffusedStat_m, rediffusedStat_m, OpAddAssign());
- reduce(stoppedPartStat_m, stoppedPartStat_m, OpAddAssign());
- if (locPartsInMat_m + bunchToMatStat_m + rediffusedStat_m + stoppedPartStat_m > 0) {
+ unsigned int tmp[4] = {locPartsInMat,
+ bunchToMatStat_m,
+ rediffusedStat_m,
+ stoppedPartStat_m};
+
+ reduce(tmp, tmp + 4, tmp, OpAddAssign());
+
+ totalPartsInMat_m = tmp[0];
+ bunchToMatStat_m = tmp[1];
+ rediffusedStat_m = tmp[2];
+ stoppedPartStat_m = tmp[3];
+
+ if (totalPartsInMat_m > 0 ||
+ bunchToMatStat_m > 0 ||
+ rediffusedStat_m > 0 ||
+ stoppedPartStat_m > 0) {
+
OPALTimer::Timer time;
msg << level2
<< "--- CollimatorPhysics - Name " << element_ref_m->getName()
@@ -685,18 +617,15 @@ void CollimatorPhysics::print(Inform &msg) {
<< std::setw(21) << "entered: " << Util::toStringWithThousandSep(bunchToMatStat_m) << "\n"
<< std::setw(21) << "rediffused: " << Util::toStringWithThousandSep(rediffusedStat_m) << "\n"
<< std::setw(21) << "stopped: " << Util::toStringWithThousandSep(stoppedPartStat_m) << "\n"
- << std::setw(21) << "total in material: " << Util::toStringWithThousandSep(locPartsInMat_m)
+ << std::setw(21) << "total in material: " << Util::toStringWithThousandSep(totalPartsInMat_m)
<< endl;
- // msg << "Coll/Deg statistics: "
- // << " bunch to material " << bunchToMatStat_m << " rediffused " << rediffusedStat_m
- // << " stopped " << stoppedPartStat_m << endl;
}
msg.flags(ff);
}
bool CollimatorPhysics::stillActive() {
- return locPartsInMat_m != 0;
+ return totalPartsInMat_m != 0;
}
bool CollimatorPhysics::stillAlive(PartBunchBase<double, 3> *bunch) {
@@ -704,7 +633,7 @@ bool CollimatorPhysics::stillAlive(PartBunchBase<double, 3> *bunch) {
bool degraderAlive = true;
//free GPU memory in case element is degrader, it is empty and bunch has moved past it
- if (collshape_m == ElementBase::DEGRADER && locPartsInMat_m == 0) {
+ if (collshape_m == ElementBase::DEGRADER && totalPartsInMat_m == 0) {
Degrader *deg = static_cast<Degrader *>(element_ref_m);
//get the size of the degrader
@@ -759,6 +688,54 @@ void CollimatorPhysics::deleteParticleFromLocalVector() {
}
}
+void CollimatorPhysics::push() {
+ for (size_t i = 0; i < locParts_m.size(); ++i) {
+ Vector_t &R = locParts_m[i].Rincol;
+ Vector_t &P = locParts_m[i].Pincol;
+ double gamma = Util::getGamma(P);
+
+ R += 0.5 * dT_m * Physics::c * P / gamma;
+ }
+}
+
+// adjust the time step for those particles that will rediffuse to
+// vacuum such that it corresponds to the time needed to reach the
+// end of the material. Only use this time step for the computation
+// of the interaction with the material, not for the integration of
+// the particles. This will ensure that the momenta of all particles
+// are reduced by approcimately the same amount in computeEnergyLoss.
+void CollimatorPhysics::setTimeStepForLeavingParticles() {
+ const double elementLength = element_ref_m->getElementLength();
+
+ for (size_t i = 0; i < locParts_m.size(); ++i) {
+ Vector_t &R = locParts_m[i].Rincol;
+ Vector_t &P = locParts_m[i].Pincol;
+ double &dt = locParts_m[i].DTincol;
+ double gamma = Util::getGamma(P);
+ Vector_t stepLength = dT_m * Physics::c * P / gamma;
+
+ if (R(2) < elementLength &&
+ R(2) + stepLength(2) > elementLength) {
+ // particle is likely to leave material
+ double distance = elementLength - R(2);
+ double tau = distance / stepLength(2);
+
+ PAssert_LE(tau, 1.0);
+ PAssert_GE(tau, 0.0);
+
+ dt *= tau;
+ }
+ }
+}
+
+void CollimatorPhysics::resetTimeStep() {
+ for (size_t i = 0; i < locParts_m.size(); ++i) {
+ double &dt = locParts_m[i].DTincol;
+ dt = dT_m;
+ }
+
+}
+
#ifdef OPAL_DKS
namespace {
@@ -767,6 +744,94 @@ namespace {
}
}
+void CollimatorPhysics::applyDKS(PartBunchBase<double, 3> *bunch,
+ const std::pair<Vector_t, double> &boundingSphere,
+ size_t numParticlesInSimulation) {
+ bool onlyOneLoopOverParticles = ! (allParticleInMat_m);
+
+ //if firs call to apply setup needed accelerator resources
+ setupCollimatorDKS(bunch, numParticlesInSimulation);
+
+ do {
+ IpplTimings::startTimer(DegraderLoopTimer_m);
+
+ //write particles to GPU if there are any to write
+ if (dksParts_m.size() > 0) {
+ //wrtie data from dksParts_m to the end of mem_ptr (offset = numparticles)
+ dksbase.writeDataAsync<PART_DKS>(mem_ptr, &dksParts_m[0],
+ dksParts_m.size(), -1, numparticles);
+
+ //update number of particles on Device
+ numparticles += dksParts_m.size();
+
+ //free locParts_m vector
+ dksParts_m.erase(dksParts_m.begin(), dksParts_m.end());
+ }
+
+ //execute CollimatorPhysics kernels on GPU if any particles are there
+ if (numparticles > 0) {
+ dksbase.callCollimatorPhysics2(mem_ptr, par_ptr, numparticles);
+ }
+
+ //sort device particles and get number of particles comming back to bunch
+ int numaddback = 0;
+ if (numparticles > 0) {
+ dksbase.callCollimatorPhysicsSort(mem_ptr, numparticles, numaddback);
+ }
+
+ //read particles from GPU if any are comming out of material
+ if (numaddback > 0) {
+
+ //resize dksParts_m to hold particles that need to go back to bunch
+ dksParts_m.resize(numaddback);
+
+ //read particles that need to be added back to bunch
+ //particles that need to be added back are at the end of Device array
+ dksbase.readData<PART_DKS>(mem_ptr, &dksParts_m[0], numaddback,
+ numparticles - numaddback);
+
+ //add particles back to the bunch
+ for (unsigned int i = 0; i < dksParts_m.size(); ++i) {
+ if (dksParts_m[i].label == -2) {
+ addBackToBunchDKS(bunch, i);
+ rediffusedStat_m++;
+ } else {
+ stoppedPartStat_m++;
+ lossDs_m->addParticle(dksParts_m[i].Rincol, dksParts_m[i].Pincol,
+ -locParts_m[dksParts_m[i].localID].IDincol);
+ }
+ }
+
+ //erase particles that came from device from host array
+ dksParts_m.erase(dksParts_m.begin(), dksParts_m.end());
+
+ //update number of particles on Device
+ numparticles -= numaddback;
+ }
+
+ IpplTimings::stopTimer(DegraderLoopTimer_m);
+
+ if (onlyOneLoopOverParticles)
+ copyFromBunchDKS(bunch, boundingSphere);
+
+ T_m += dT_m;
+
+ totalPartsInMat_m = numparticles;
+ reduce(totalPartsInMat_m, totalPartsInMat_m, OpAddAssign());
+
+ int maxPerNode = bunch->getLocalNum();
+ reduce(maxPerNode, maxPerNode, OpMaxAssign());
+
+ //more than one loop only if all the particles are in this degrader
+ if (allParticleInMat_m) {
+ onlyOneLoopOverParticles = (totalPartsInMat_m <= 0);
+ } else {
+ onlyOneLoopOverParticles = true;
+ }
+
+ } while (onlyOneLoopOverParticles == false);
+}
+
void CollimatorPhysics::addBackToBunchDKS(PartBunchBase<double, 3> *bunch, unsigned i) {
int id = dksParts_m[i].localID;
@@ -925,4 +990,4 @@ void CollimatorPhysics::deleteParticleFromLocalVectorDKS() {
}
-#endif
+#endif
\ No newline at end of file
diff --git a/src/Classic/Solvers/CollimatorPhysics.hh b/src/Classic/Solvers/CollimatorPhysics.hh
index 87ac5b8c7..169c2a76d 100644
--- a/src/Classic/Solvers/CollimatorPhysics.hh
+++ b/src/Classic/Solvers/CollimatorPhysics.hh
@@ -78,7 +78,10 @@ struct InsideTester {
class CollimatorPhysics: public ParticleMatterInteractionHandler {
public:
- CollimatorPhysics(const std::string &name, ElementBase *element, std::string &mat);
+ CollimatorPhysics(const std::string &name,
+ ElementBase *element,
+ std::string &mat,
+ bool enableRutherford);
~CollimatorPhysics();
void apply(PartBunchBase<double, 3> *bunch,
@@ -96,28 +99,28 @@ public:
size_t getParticlesInMat();
unsigned getRediffused();
unsigned int getNumEntered();
- void doPhysics(PartBunchBase<double, 3> *bunch);
+ void computeInteraction();
- virtual bool computeEnergyLoss(double &Eng, const double deltat, bool includeFluctuations = true) const;
+ virtual bool computeEnergyLoss(Vector_t &P,
+ const double deltat,
+ bool includeFluctuations = true) const;
private:
void configureMaterialParameters();
- void computeCoulombScattering(Vector_t &R, Vector_t &P, const double &deltat);
+ void computeCoulombScattering(Vector_t &R,
+ Vector_t &P,
+ double dt);
- void applyRotation(double &px,
- double &pz,
- double &x,
- double &z,
+ void applyRotation(Vector_t &P,
+ Vector_t &R,
double xplane,
- double Norm_P,
- double thetacou,
- double deltas,
- int coord);
+ double thetacou);
+ void applyRandomRotation(Vector_t &P, double theta0);
void copyFromBunch(PartBunchBase<double, 3> *bunch,
const std::pair<Vector_t, double> &boundingSphere);
- void addBackToBunch(PartBunchBase<double, 3> *bunch, unsigned i);
+ void addBackToBunch(PartBunchBase<double, 3> *bunch);
void applyNonDKS(PartBunchBase<double, 3> *bunch,
const std::pair<Vector_t, double> &boundingSphere,
@@ -128,7 +131,7 @@ private:
size_t numParticlesInSimulation);
void copyFromBunchDKS(PartBunchBase<double, 3> *bunch,
const std::pair<Vector_t, double> &boundingSphere);
- void addBackToBunchDKS(PartBunchBase<double, 3> *bunch, unsigned i);
+ void addBackToBunchDKS(PartBunchBase<double, 3> *bunch);
void setupCollimatorDKS(PartBunchBase<double, 3> *bunch, size_t numParticlesInSimulation);
void clearCollimatorDKS();
@@ -141,6 +144,10 @@ private:
void calcStat(double Eng);
+ void push();
+ void resetTimeStep();
+ void setTimeStepForLeavingParticles();
+
double T_m; // own time, maybe larger than in the bunch object
double dT_m; // dt from bunch
@@ -169,8 +176,8 @@ private:
unsigned int stoppedPartStat_m;
// number of particles that leave the material in current step
unsigned int rediffusedStat_m;
- // number of local particles that are in the material
- size_t locPartsInMat_m;
+ // total number of particles that are in the material
+ unsigned int totalPartsInMat_m;
// some statistics
double Eavg_m; // average kinetic energy
@@ -181,6 +188,7 @@ private:
std::unique_ptr<LossDataSink> lossDs_m;
+ bool enableRutherford_m;
#ifdef OPAL_DKS
DKSOPAL dksbase;
int curandInitSet;
@@ -218,7 +226,7 @@ double CollimatorPhysics::getTime() {
inline
size_t CollimatorPhysics::getParticlesInMat() {
- return locPartsInMat_m;
+ return totalPartsInMat_m;
}
inline
diff --git a/src/Classic/Solvers/ParticleMatterInteractionHandler.hh b/src/Classic/Solvers/ParticleMatterInteractionHandler.hh
index f928a8755..de4127968 100644
--- a/src/Classic/Solvers/ParticleMatterInteractionHandler.hh
+++ b/src/Classic/Solvers/ParticleMatterInteractionHandler.hh
@@ -31,7 +31,9 @@ public:
void updateElement(ElementBase *newref);
ElementBase* getElement();
- virtual bool computeEnergyLoss(double &Eng, const double deltat, bool includeFluctuations = true) const = 0;
+ virtual bool computeEnergyLoss(Vector_t &P,//double &Eng,
+ const double deltat,
+ bool includeFluctuations = true) const = 0;
protected:
ElementBase *element_ref_m;
bool allParticleInMat_m; ///< if all particles are in matter stay inside the particle matter interaction
diff --git a/src/Structure/ParticleMatterInteraction.cpp b/src/Structure/ParticleMatterInteraction.cpp
index 8acb90b41..74602edc5 100644
--- a/src/Structure/ParticleMatterInteraction.cpp
+++ b/src/Structure/ParticleMatterInteraction.cpp
@@ -37,10 +37,7 @@ namespace {
// DESCRIPTION OF SINGLE PARTICLE:
TYPE, // The type of the wake
MATERIAL, // From of the tube
- RADIUS, // Radius of the tube
- SIGMA,
- TAU,
- NPART,
+ ENABLERUTHERFORD,
SIZE
};
}
@@ -56,16 +53,8 @@ ParticleMatterInteraction::ParticleMatterInteraction():
itsAttr[MATERIAL] = Attributes::makeString
("MATERIAL", "The material of the surface");
- itsAttr[RADIUS] = Attributes::makeReal
- ("RADIUS", "The radius of the beam pipe [m]");
-
- itsAttr[SIGMA] = Attributes::makeReal
- ("SIGMA", "Material constant dependant on the beam pipe material");
-
- itsAttr[TAU] = Attributes::makeReal
- ("TAU", "Material constant dependant on the beam pipe material");
-
- itsAttr[NPART] = Attributes::makeReal("NPART", "Number of particles in bunch");
+ itsAttr[ENABLERUTHERFORD] = Attributes::makeBool("ENABLERUTHERFORD",
+ "Enable large angle scattering", true);
ParticleMatterInteraction *defParticleMatterInteraction = clone("UNNAMED_PARTICLEMATTERINTERACTION");
defParticleMatterInteraction->builtin = true;
@@ -130,15 +119,15 @@ void ParticleMatterInteraction::initParticleMatterInteractionHandler(ElementBase
*gmsg << "* ParticleMatterInteraction::initParticleMatterInteractionHandler " << endl;
*gmsg << "* ********************************************************************************** " << endl;
- itsElement_m = &element;
- material_m = Util::toUpper(Attributes::getString(itsAttr[MATERIAL]));
+ std::string material = Util::toUpper(Attributes::getString(itsAttr[MATERIAL]));
+ bool enableRutherford = Attributes::getBool(itsAttr[ENABLERUTHERFORD]);
const std::string type = Util::toUpper(Attributes::getString(itsAttr[TYPE]));
if(type == "CCOLLIMATOR" ||
type == "COLLIMATOR" ||
type == "DEGRADER") {
- handler_m = new CollimatorPhysics(getOpalName(), itsElement_m, material_m);
+ handler_m = new CollimatorPhysics(getOpalName(), &element, material, enableRutherford);
*gmsg << *this << endl;
} else {
handler_m = 0;
diff --git a/src/Structure/ParticleMatterInteraction.h b/src/Structure/ParticleMatterInteraction.h
index 491caf254..84b97ce40 100644
--- a/src/Structure/ParticleMatterInteraction.h
+++ b/src/Structure/ParticleMatterInteraction.h
@@ -72,17 +72,6 @@ private:
// Clone constructor.
ParticleMatterInteraction(const std::string &name, ParticleMatterInteraction *parent);
-
- // The particle reference data.
- PartData reference;
-
- // The conversion from GeV to eV.
- static const double energy_scale;
-
- // the element the particle mater interaction is attached to
- ElementBase *itsElement_m;
- std::string material_m;
-
};
inline std::ostream &operator<<(std::ostream &os, const ParticleMatterInteraction &b) {
@@ -90,4 +79,4 @@ inline std::ostream &operator<<(std::ostream &os, const ParticleMatterInteractio
return os;
}
-#endif // OPAL_PARTICLEMATTERINTERACTION_HH
+#endif // OPAL_PARTICLEMATTERINTERACTION_HH
\ No newline at end of file
diff --git a/tools/emacs/opal.el b/tools/emacs/opal.el
index 31a41c02e..c2efbfaf4 100644
--- a/tools/emacs/opal.el
+++ b/tools/emacs/opal.el
@@ -138,11 +138,11 @@
)
"Highlighting expressions for OPAL mode (seqediting).")
-;(concat "\\<" (regexp-opt '("A" "ADD" "AMPLITUDE_MODEL" "ANGLE" "APERTURE" "APVETO" "AT" "AUTOPHASE" "B" "BBOXINCR" "BCFFTZ" "BCFFTT" "BCFFTX" "BCFFTY" "BCURRENT" "BEAM_PHIINIT" "BEAM_PRINIT" "BEAM_RINIT" "BFREQ" "BIRTH_CONTROL" "BMAX" "BOUNDPDESTROYFQ" "BSCALE" "BY" "CALLS" "CATHTEMP" "CENTRE" "CHARGE" "CLASS" "CLEAR" "CMD" "COEFDENOM" "COEFDENOMPHI" "COEFNUM" "COEFNUMPHI" "COLUMN" "CONDUCT" "CONSTRAINTS" "CONST_LENGTH" "CONV_HVOL_PROG" "CORRX" "CORRY" "CORRZ" "CSRDUMP" "CROSSOVER" "CUTOFFLONG" "CUTOFFPX" "CUTOFFPY" "CUTOFFPZ" "CUTOFFR" "CUTOFFX" "CUTOFFY" "CYHARMON" "CZERO" "DESIGNENERGY" "DK1" "DK1S" "DK2" "DK2S" "DKN" "DKNR" "DKS" "DKSR" "DLAG" "DPHI" "DPSI" "DS" "DT" "DTHETA" "DUMP" "DUMP_DAT" "DUMP_FREQ" "DUMP_OFFSPRING" "DVARS" "DVOLT" "DX" "DY" "DZ" "E1" "E2" "EBDUMP" "ECHO" "EKIN" "ELASER" "ELEMEDGE" "ELEMENT" "EMISSIONMODEL" "EMISSIONSTEPS" "EMITTED" "ENABLEHDF5" "ENERGY" "ENDSEQUENCE" "END_NORMAL_X" "END_NORMAL_Y" "END_POSITION_X" "END_POSITION_Y" "EPSILON" "ESCALE" "ET" "EVERYSTEP" "EX" "EXPECTED_HYPERVOL" "EXPR" "EY" "FE" "FGEOM" "FIELDMAPDIR" "FIELDMAPDIR" "FILE" "FINT" "FMAPFN" "FNAME" "FORM" "FREQ" "FREQUENCY_MODEL" "FROM" "FSTYPE" "FTOSCAMPLITUDE" "FTOSCPERIODS" "FULL" "GAMMA" "GAP" "GAPWIDTH" "GENE_MUTATION_PROBABILITY" "GEOMETRY" "GREENSF" "H1" "H2" "HAPERT" "HARMON" "HARMONIC_NUMBER" "HGAP" "HKICK" "HYPERVOLREFERENCE" "IDEALIZED" "INITIAL_OPTIMIZATION" "IMAGENAME" "INFO" "INITIALPOPULATION" "INPUT" "INPUTMOUNITS" "INTENSITYCUT" "INTERPL" "IS_CLOSED" "ITSOLVER" "K0" "K0S" "K1" "K1S" "K2" "K2S" "K3" "K3S" "KEYWORD" "KICK" "KN" "KS" "L" "LAG" "LASERPROFFN" "LATTICE_PHIINIT" "LATTICE_RINIT" "LATTICE_THETAINIT" "LENGTH" "LEVEL" "LOGBENDTRAJECTORY" "LOWER" "LOWERBOUND" "MASS" "MAXGENERATIONS" "MAXITERS" "MAXR" "MAXSTEPS" "MAXZ" "METHOD" "MINR" "MINZ" "MODE" "MREX" "MREY" "MSCALX" "MSCALY" "MT" "MUTATION_PROBABILITY" "MX" "MY" "NBIN" "NFREQ" "NLEFT" "NO" "NPART" "NPOINTS" "NRIGHT" "NUMCELLS" "NUM_COWORKERS" "NUM_IND_GEN" "NUM_MASTERS" "OBJECTIVES" "OFFSETPX" "OFFSETPY" "OFFSETPZ" "OFFSETT" "OFFSETX" "OFFSETY" "OFFSETZ" "ONE_PILOT_CONVERGE" "OPCHARGE" "OPMASS" "OPYIELD" "ORDER" "ORIENTATION" "ORIGIN" "OUTDIR" "OUTFN" "OUTPUT" "P1" "P2" "P3" "P4" "PARFFTT" "PARFFTX" "PARFFTY" "PARTICLE" "PARTICLEMATTERINTERACTION" "PATTERN" "PC" "PDIS" "PHI" "PHI0" "PHIINIT" "PHIMAX" "PHIMIN" "PLANE" "POLYORDER" "PRECMODE" "PRINIT" "PSDUMPEACHTURN" "PSDUMPFRAME" "PSDUMPFREQ" "PSI" "PTC" "PYMULT" "PZINIT" "PZMULT" "R51" "R52" "R61" "R62" "RADIUS" "RANDOM" "RANGE" "RASTER" "RECOMBINATION_PROBABILITY" "REFER" "REFPOS" "REPARTFREQ" "RESET" "RFFREQ" "RFMAPFN" "RFPHI" "RINIT" "RMAX" "RMIN" "ROTATION" "ROW" "S" "SAMPLINGS" "SCALABLE" "SEED" "SELECTED" "SEQUENCE" "SIGMA" "SIGMAPX" "SIGMAPY" "SIGMAPZ" "SIGMAR" "SIGMAT" "SIGMAX" "SIGMAY" "SIGMAZ" "SIGX" "SIGY" "SIMBIN_CROSSOVER_NU" "SIMTMPDIR" "SLPTC" "SOL_SYNCH" "SPLIT" "SPTDUMPFREQ" "STARTPOPULATION" "STATDUMPFREQ" "STEP" "STEPSPERTURN" "STOP" "STRING" "SUPERPOSE" "SYMMETRY" "T0" "TABLE" "TAU" "TELL" "TEMPLATEDIR" "TFALL" "THETA" "THIN" "THRESHOLD" "TIME" "TIMEINTEGRATOR" "TMULT" "TO" "TOL" "TOLERANCE" "TPULSEFWHM" "TRACE" "TRIMCOILTHRESHOLD" "TRISE" "TURNS" "TYPE" "UPPER" "UPPERBOUND" "VARIABLE" "VERIFY" "VERSION" "VKICK" "VMAX" "VMIN" "VOLT" "W" "WAKEF" "WARN" "WARP" "WEIGHT" "WIDTH" "WRITETOFILE" "X" "XEND" "XMA" "XMULT" "XSIZE" "XSTART" "Y" "YEND" "YMA" "YMULT" "YSIZE" "YSTART" "Z" "Z0" "ZEND" "ZINIT" "ZSTART" "ZSTOP") t) "\\>")
+;(concat "\\<" (regexp-opt '("A" "ADD" "AMPLITUDE_MODEL" "ANGLE" "APERTURE" "APVETO" "AT" "AUTOPHASE" "B" "BBOXINCR" "BCFFTZ" "BCFFTT" "BCFFTX" "BCFFTY" "BCURRENT" "BEAM_PHIINIT" "BEAM_PRINIT" "BEAM_RINIT" "BFREQ" "BIRTH_CONTROL" "BMAX" "BOUNDPDESTROYFQ" "BSCALE" "BY" "CALLS" "CATHTEMP" "CENTRE" "CHARGE" "CLASS" "CLEAR" "CMD" "COEFDENOM" "COEFDENOMPHI" "COEFNUM" "COEFNUMPHI" "COLUMN" "CONDUCT" "CONSTRAINTS" "CONST_LENGTH" "CONV_HVOL_PROG" "CORRX" "CORRY" "CORRZ" "CSRDUMP" "CROSSOVER" "CUTOFFLONG" "CUTOFFPX" "CUTOFFPY" "CUTOFFPZ" "CUTOFFR" "CUTOFFX" "CUTOFFY" "CYHARMON" "CZERO" "DESIGNENERGY" "DK1" "DK1S" "DK2" "DK2S" "DKN" "DKNR" "DKS" "DKSR" "DLAG" "DPHI" "DPSI" "DS" "DT" "DTHETA" "DUMP" "DUMP_DAT" "DUMP_FREQ" "DUMP_OFFSPRING" "DVARS" "DVOLT" "DX" "DY" "DZ" "E1" "E2" "EBDUMP" "ECHO" "EKIN" "ELASER" "ELEMEDGE" "ELEMENT" "EMISSIONMODEL" "EMISSIONSTEPS" "EMITTED" "ENABLEHDF5" "ENABLERUTHERFORD" "ENERGY" "ENDSEQUENCE" "END_NORMAL_X" "END_NORMAL_Y" "END_POSITION_X" "END_POSITION_Y" "EPSILON" "ESCALE" "ET" "EVERYSTEP" "EX" "EXPECTED_HYPERVOL" "EXPR" "EY" "FE" "FGEOM" "FIELDMAPDIR" "FIELDMAPDIR" "FILE" "FINT" "FMAPFN" "FNAME" "FORM" "FREQ" "FREQUENCY_MODEL" "FROM" "FSTYPE" "FTOSCAMPLITUDE" "FTOSCPERIODS" "FULL" "GAMMA" "GAP" "GAPWIDTH" "GENE_MUTATION_PROBABILITY" "GEOMETRY" "GREENSF" "H1" "H2" "HAPERT" "HARMON" "HARMONIC_NUMBER" "HGAP" "HKICK" "HYPERVOLREFERENCE" "IDEALIZED" "INITIAL_OPTIMIZATION" "IMAGENAME" "INFO" "INITIALPOPULATION" "INPUT" "INPUTMOUNITS" "INTENSITYCUT" "INTERPL" "IS_CLOSED" "ITSOLVER" "K0" "K0S" "K1" "K1S" "K2" "K2S" "K3" "K3S" "KEYWORD" "KICK" "KN" "KS" "L" "LAG" "LASERPROFFN" "LATTICE_PHIINIT" "LATTICE_RINIT" "LATTICE_THETAINIT" "LENGTH" "LEVEL" "LOGBENDTRAJECTORY" "LOWER" "LOWERBOUND" "MASS" "MATERIAL" "MAXGENERATIONS" "MAXITERS" "MAXR" "MAXSTEPS" "MAXZ" "METHOD" "MINR" "MINZ" "MODE" "MREX" "MREY" "MSCALX" "MSCALY" "MT" "MUTATION_PROBABILITY" "MX" "MY" "NBIN" "NFREQ" "NLEFT" "NO" "NPART" "NPOINTS" "NRIGHT" "NUMCELLS" "NUM_COWORKERS" "NUM_IND_GEN" "NUM_MASTERS" "OBJECTIVES" "OFFSETPX" "OFFSETPY" "OFFSETPZ" "OFFSETT" "OFFSETX" "OFFSETY" "OFFSETZ" "ONE_PILOT_CONVERGE" "OPCHARGE" "OPMASS" "OPYIELD" "ORDER" "ORIENTATION" "ORIGIN" "OUTDIR" "OUTFN" "OUTPUT" "P1" "P2" "P3" "P4" "PARFFTT" "PARFFTX" "PARFFTY" "PARTICLE" "PARTICLEMATTERINTERACTION" "PATTERN" "PC" "PDIS" "PHI" "PHI0" "PHIINIT" "PHIMAX" "PHIMIN" "PLANE" "POLYORDER" "PRECMODE" "PRINIT" "PSDUMPEACHTURN" "PSDUMPFRAME" "PSDUMPFREQ" "PSI" "PTC" "PYMULT" "PZINIT" "PZMULT" "R51" "R52" "R61" "R62" "RADIUS" "RANDOM" "RANGE" "RASTER" "RECOMBINATION_PROBABILITY" "REFER" "REFPOS" "REPARTFREQ" "RESET" "RFFREQ" "RFMAPFN" "RFPHI" "RINIT" "RMAX" "RMIN" "ROTATION" "ROW" "S" "SAMPLINGS" "SCALABLE" "SEED" "SELECTED" "SEQUENCE" "SIGMA" "SIGMAPX" "SIGMAPY" "SIGMAPZ" "SIGMAR" "SIGMAT" "SIGMAX" "SIGMAY" "SIGMAZ" "SIGX" "SIGY" "SIMBIN_CROSSOVER_NU" "SIMTMPDIR" "SLPTC" "SOL_SYNCH" "SPLIT" "SPTDUMPFREQ" "STARTPOPULATION" "STATDUMPFREQ" "STEP" "STEPSPERTURN" "STOP" "STRING" "SUPERPOSE" "SYMMETRY" "T0" "TABLE" "TAU" "TELL" "TEMPLATEDIR" "TFALL" "THETA" "THIN" "THRESHOLD" "TIME" "TIMEINTEGRATOR" "TMULT" "TO" "TOL" "TOLERANCE" "TPULSEFWHM" "TRACE" "TRIMCOILTHRESHOLD" "TRISE" "TURNS" "TYPE" "UPPER" "UPPERBOUND" "VARIABLE" "VERIFY" "VERSION" "VKICK" "VMAX" "VMIN" "VOLT" "W" "WAKEF" "WARN" "WARP" "WEIGHT" "WIDTH" "WRITETOFILE" "X" "XEND" "XMA" "XMULT" "XSIZE" "XSTART" "Y" "YEND" "YMA" "YMULT" "YSIZE" "YSTART" "Z" "Z0" "ZEND" "ZINIT" "ZSTART" "ZSTOP") t) "\\>")
(defconst opal-font-lock-keywords-parameters
(list
- '("\\<\\(A\\(?:DD\\|MPLITUDE_MODEL\\|NGLE\\|P\\(?:ERTURE\\|VETO\\)\\|T\\|UTOPHASE\\)\\|B\\(?:BOXINCR\\|C\\(?:FFT[TXYZ]\\|URRENT\\)\\|EAM_\\(?:\\(?:P\\(?:HI\\|R\\)\\|R\\)INIT\\)\\|FREQ\\|IRTH_CONTROL\\|MAX\\|OUNDPDESTROYFQ\\|SCALE\\|Y\\)\\|C\\(?:A\\(?:LLS\\|THTEMP\\)\\|ENTRE\\|HARGE\\|L\\(?:ASS\\|EAR\\)\\|MD\\|O\\(?:EF\\(?:DENOM\\(?:PHI\\)?\\|NUM\\(?:PHI\\)?\\)\\|LUMN\\|N\\(?:DUCT\\|ST\\(?:RAINTS\\|_LENGTH\\)\\|V_HVOL_PROG\\)\\|RR[XYZ]\\)\\|ROSSOVER\\|SRDUMP\\|UTOFF\\(?:LONG\\|P[XYZ]\\|[RXY]\\)\\|YHARMON\\|ZERO\\)\\|D\\(?:ESIGNENERGY\\|K\\(?:1S\\|2S\\|[NS]R\\|[12NS]\\)\\|LAG\\|P\\(?:[HS]I\\)\\|THETA\\|UMP\\(?:_\\(?:DAT\\|FREQ\\|OFFSPRING\\)\\)?\\|V\\(?:ARS\\|OLT\\)\\|[STXYZ]\\)\\|E\\(?:BDUMP\\|CHO\\|KIN\\|L\\(?:ASER\\|EME\\(?:DGE\\|NT\\)\\)\\|MI\\(?:SSION\\(?:MODEL\\|STEPS\\)\\|TTED\\)\\|N\\(?:ABLEHDF5\\|D\\(?:SEQUENCE\\|_\\(?:NORMAL_[XY]\\|POSITION_[XY]\\)\\)\\|ERGY\\)\\|PSILON\\|SCALE\\|VERYSTEP\\|XP\\(?:ECTED_HYPERVOL\\|R\\)\\|[12TXY]\\)\\|F\\(?:E\\|GEOM\\|I\\(?:ELDMAPDIR\\|LE\\|NT\\)\\|MAPFN\\|NAME\\|ORM\\|R\\(?:EQ\\(?:UENCY_MODEL\\)?\\|OM\\)\\|STYPE\\|TOSC\\(?:AMPLITUDE\\|PERIODS\\)\\|ULL\\)\\|G\\(?:A\\(?:MMA\\|P\\(?:WIDTH\\)?\\)\\|E\\(?:\\(?:NE_MUTATION_PROBABILIT\\|OMETR\\)Y\\)\\|REENSF\\)\\|H\\(?:A\\(?:PERT\\|RMON\\(?:IC_NUMBER\\)?\\)\\|GAP\\|KICK\\|YPERVOLREFERENCE\\|[12]\\)\\|I\\(?:DEALIZED\\|MAGENAME\\|N\\(?:FO\\|ITIAL\\(?:\\(?:POPUL\\|_OPTIMIZ\\)ATION\\)\\|PUT\\(?:MOUNITS\\)?\\|TE\\(?:NSITYCUT\\|RPL\\)\\)\\|S_CLOSED\\|TSOLVER\\)\\|K\\(?:0S\\|1S\\|2S\\|3S\\|EYWORD\\|ICK\\|[0-3NS]\\)\\|L\\(?:A\\(?:G\\|SERPROFFN\\|TTICE_\\(?:\\(?:PHI\\|R\\|THETA\\)INIT\\)\\)\\|E\\(?:NGTH\\|VEL\\)\\|O\\(?:GBENDTRAJECTORY\\|WER\\(?:BOUND\\)?\\)\\)\\|M\\(?:A\\(?:SS\\|X\\(?:GENERATIONS\\|ITERS\\|STEPS\\|[RZ]\\)\\)\\|ETHOD\\|IN[RZ]\\|ODE\\|RE[XY]\\|SCAL[XY]\\|UTATION_PROBABILITY\\|[TXY]\\)\\|N\\(?:BIN\\|FREQ\\|LEFT\\|O\\|P\\(?:ART\\|OINTS\\)\\|RIGHT\\|UM\\(?:CELLS\\|_\\(?:COWORKERS\\|IND_GEN\\|MASTERS\\)\\)\\)\\|O\\(?:BJECTIVES\\|FFSET\\(?:P[XYZ]\\|[TXYZ]\\)\\|NE_PILOT_CONVERGE\\|P\\(?:CHARGE\\|MASS\\|YIELD\\)\\|R\\(?:DER\\|I\\(?:\\(?:ENTATIO\\|GI\\)N\\)\\)\\|UT\\(?:DIR\\|FN\\|PUT\\)\\)\\|P\\(?:A\\(?:R\\(?:FFT[TXY]\\|TICLE\\(?:MATTERINTERACTION\\)?\\)\\|TTERN\\)\\|DIS\\|HI\\(?:0\\|INIT\\|M\\(?:AX\\|IN\\)\\)?\\|LANE\\|OLYORDER\\|R\\(?:ECMODE\\|INIT\\)\\|S\\(?:DUMP\\(?:EACHTURN\\|FR\\(?:AME\\|EQ\\)\\)\\|I\\)\\|TC\\|\\(?:YMUL\\|Z\\(?:INI\\|MUL\\)\\)T\\|[1-4C]\\)\\|R\\(?:5[12]\\|6[12]\\|A\\(?:DIUS\\|N\\(?:DOM\\|GE\\)\\|STER\\)\\|E\\(?:COMBINATION_PROBABILITY\\|F\\(?:ER\\|POS\\)\\|PARTFREQ\\|SET\\)\\|F\\(?:FREQ\\|MAPFN\\|PHI\\)\\|INIT\\|M\\(?:AX\\|IN\\)\\|O\\(?:TATION\\|W\\)\\)\\|S\\(?:AMPLINGS\\|CALABLE\\|E\\(?:ED\\|LECTED\\|QUENCE\\)\\|I\\(?:G\\(?:MA\\(?:P[XYZ]\\|[RTXYZ]\\)?\\|[XY]\\)\\|M\\(?:BIN_CROSSOVER_NU\\|TMPDIR\\)\\)\\|LPTC\\|OL_SYNCH\\|P\\(?:LIT\\|TDUMPFREQ\\)\\|T\\(?:A\\(?:RTPOPULATION\\|TDUMPFREQ\\)\\|EP\\(?:SPERTURN\\)?\\|OP\\|RING\\)\\|UPERPOSE\\|YMMETRY\\)\\|T\\(?:A\\(?:BLE\\|U\\)\\|E\\(?:LL\\|MPLATEDIR\\)\\|FALL\\|H\\(?:ETA\\|IN\\|RESHOLD\\)\\|IME\\(?:INTEGRATOR\\)?\\|MULT\\|OL\\(?:ERANCE\\)?\\|PULSEFWHM\\|R\\(?:ACE\\|I\\(?:MCOILTHRESHOLD\\|SE\\)\\)\\|URNS\\|YPE\\|[0O]\\)\\|UPPER\\(?:BOUND\\)?\\|V\\(?:ARIABLE\\|ER\\(?:IFY\\|SION\\)\\|KICK\\|M\\(?:AX\\|IN\\)\\|OLT\\)\\|W\\(?:A\\(?:KEF\\|R[NP]\\)\\|EIGHT\\|IDTH\\|RITETOFILE\\)\\|X\\(?:END\\|M\\(?:A\\|ULT\\)\\|S\\(?:IZE\\|TART\\)\\)\\|Y\\(?:END\\|M\\(?:A\\|ULT\\)\\|S\\(?:IZE\\|TART\\)\\)\\|Z\\(?:0\\|END\\|INIT\\|ST\\(?:ART\\|OP\\)\\)\\|[ABLSW-Z]\\)\\>"
+ '("\\<\\(A\\(?:DD\\|MPLITUDE_MODEL\\|NGLE\\|P\\(?:ERTURE\\|VETO\\)\\|T\\|UTOPHASE\\)\\|B\\(?:BOXINCR\\|C\\(?:FFT[TXYZ]\\|URRENT\\)\\|EAM_\\(?:\\(?:P\\(?:HI\\|R\\)\\|R\\)INIT\\)\\|FREQ\\|IRTH_CONTROL\\|MAX\\|OUNDPDESTROYFQ\\|SCALE\\|Y\\)\\|C\\(?:A\\(?:LLS\\|THTEMP\\)\\|ENTRE\\|HARGE\\|L\\(?:ASS\\|EAR\\)\\|MD\\|O\\(?:EF\\(?:DENOM\\(?:PHI\\)?\\|NUM\\(?:PHI\\)?\\)\\|LUMN\\|N\\(?:DUCT\\|ST\\(?:RAINTS\\|_LENGTH\\)\\|V_HVOL_PROG\\)\\|RR[XYZ]\\)\\|ROSSOVER\\|SRDUMP\\|UTOFF\\(?:LONG\\|P[XYZ]\\|[RXY]\\)\\|YHARMON\\|ZERO\\)\\|D\\(?:ESIGNENERGY\\|K\\(?:1S\\|2S\\|[NS]R\\|[12NS]\\)\\|LAG\\|P\\(?:[HS]I\\)\\|THETA\\|UMP\\(?:_\\(?:DAT\\|FREQ\\|OFFSPRING\\)\\)?\\|V\\(?:ARS\\|OLT\\)\\|[STXYZ]\\)\\|E\\(?:BDUMP\\|CHO\\|KIN\\|L\\(?:ASER\\|EME\\(?:DGE\\|NT\\)\\)\\|MI\\(?:SSION\\(?:MODEL\\|STEPS\\)\\|TTED\\)\\|N\\(?:ABLE\\(?:HDF5\\|RUTHERFORD\\)\\|D\\(?:SEQUENCE\\|_\\(?:NORMAL_[XY]\\|POSITION_[XY]\\)\\)\\|ERGY\\)\\|PSILON\\|SCALE\\|VERYSTEP\\|XP\\(?:ECTED_HYPERVOL\\|R\\)\\|[12TXY]\\)\\|F\\(?:E\\|GEOM\\|I\\(?:ELDMAPDIR\\|LE\\|NT\\)\\|MAPFN\\|NAME\\|ORM\\|R\\(?:EQ\\(?:UENCY_MODEL\\)?\\|OM\\)\\|STYPE\\|TOSC\\(?:AMPLITUDE\\|PERIODS\\)\\|ULL\\)\\|G\\(?:A\\(?:MMA\\|P\\(?:WIDTH\\)?\\)\\|E\\(?:\\(?:NE_MUTATION_PROBABILIT\\|OMETR\\)Y\\)\\|REENSF\\)\\|H\\(?:A\\(?:PERT\\|RMON\\(?:IC_NUMBER\\)?\\)\\|GAP\\|KICK\\|YPERVOLREFERENCE\\|[12]\\)\\|I\\(?:DEALIZED\\|MAGENAME\\|N\\(?:FO\\|ITIAL\\(?:\\(?:POPUL\\|_OPTIMIZ\\)ATION\\)\\|PUT\\(?:MOUNITS\\)?\\|TE\\(?:NSITYCUT\\|RPL\\)\\)\\|S_CLOSED\\|TSOLVER\\)\\|K\\(?:0S\\|1S\\|2S\\|3S\\|EYWORD\\|ICK\\|[0-3NS]\\)\\|L\\(?:A\\(?:G\\|SERPROFFN\\|TTICE_\\(?:\\(?:PHI\\|R\\|THETA\\)INIT\\)\\)\\|E\\(?:NGTH\\|VEL\\)\\|O\\(?:GBENDTRAJECTORY\\|WER\\(?:BOUND\\)?\\)\\)\\|M\\(?:A\\(?:SS\\|TERIAL\\|X\\(?:GENERATIONS\\|ITERS\\|STEPS\\|[RZ]\\)\\)\\|ETHOD\\|IN[RZ]\\|ODE\\|RE[XY]\\|SCAL[XY]\\|UTATION_PROBABILITY\\|[TXY]\\)\\|N\\(?:BIN\\|FREQ\\|LEFT\\|O\\|P\\(?:ART\\|OINTS\\)\\|RIGHT\\|UM\\(?:CELLS\\|_\\(?:COWORKERS\\|IND_GEN\\|MASTERS\\)\\)\\)\\|O\\(?:BJECTIVES\\|FFSET\\(?:P[XYZ]\\|[TXYZ]\\)\\|NE_PILOT_CONVERGE\\|P\\(?:CHARGE\\|MASS\\|YIELD\\)\\|R\\(?:DER\\|I\\(?:\\(?:ENTATIO\\|GI\\)N\\)\\)\\|UT\\(?:DIR\\|FN\\|PUT\\)\\)\\|P\\(?:A\\(?:R\\(?:FFT[TXY]\\|TICLE\\(?:MATTERINTERACTION\\)?\\)\\|TTERN\\)\\|DIS\\|HI\\(?:0\\|INIT\\|M\\(?:AX\\|IN\\)\\)?\\|LANE\\|OLYORDER\\|R\\(?:ECMODE\\|INIT\\)\\|S\\(?:DUMP\\(?:EACHTURN\\|FR\\(?:AME\\|EQ\\)\\)\\|I\\)\\|TC\\|\\(?:YMUL\\|Z\\(?:INI\\|MUL\\)\\)T\\|[1-4C]\\)\\|R\\(?:5[12]\\|6[12]\\|A\\(?:DIUS\\|N\\(?:DOM\\|GE\\)\\|STER\\)\\|E\\(?:COMBINATION_PROBABILITY\\|F\\(?:ER\\|POS\\)\\|PARTFREQ\\|SET\\)\\|F\\(?:FREQ\\|MAPFN\\|PHI\\)\\|INIT\\|M\\(?:AX\\|IN\\)\\|O\\(?:TATION\\|W\\)\\)\\|S\\(?:AMPLINGS\\|CALABLE\\|E\\(?:ED\\|LECTED\\|QUENCE\\)\\|I\\(?:G\\(?:MA\\(?:P[XYZ]\\|[RTXYZ]\\)?\\|[XY]\\)\\|M\\(?:BIN_CROSSOVER_NU\\|TMPDIR\\)\\)\\|LPTC\\|OL_SYNCH\\|P\\(?:LIT\\|TDUMPFREQ\\)\\|T\\(?:A\\(?:RTPOPULATION\\|TDUMPFREQ\\)\\|EP\\(?:SPERTURN\\)?\\|OP\\|RING\\)\\|UPERPOSE\\|YMMETRY\\)\\|T\\(?:A\\(?:BLE\\|U\\)\\|E\\(?:LL\\|MPLATEDIR\\)\\|FALL\\|H\\(?:ETA\\|IN\\|RESHOLD\\)\\|IME\\(?:INTEGRATOR\\)?\\|MULT\\|OL\\(?:ERANCE\\)?\\|PULSEFWHM\\|R\\(?:ACE\\|I\\(?:MCOILTHRESHOLD\\|SE\\)\\)\\|URNS\\|YPE\\|[0O]\\)\\|UPPER\\(?:BOUND\\)?\\|V\\(?:ARIABLE\\|ER\\(?:IFY\\|SION\\)\\|KICK\\|M\\(?:AX\\|IN\\)\\|OLT\\)\\|W\\(?:A\\(?:KEF\\|R[NP]\\)\\|EIGHT\\|IDTH\\|RITETOFILE\\)\\|X\\(?:END\\|M\\(?:A\\|ULT\\)\\|S\\(?:IZE\\|TART\\)\\)\\|Y\\(?:END\\|M\\(?:A\\|ULT\\)\\|S\\(?:IZE\\|TART\\)\\)\\|Z\\(?:0\\|END\\|INIT\\|ST\\(?:ART\\|OP\\)\\)\\|[ABLSW-Z]\\)\\>"
. font-lock-variable-name-face)
)
"Highlighting expressions for OPAL mode (parameters).")
--
GitLab