Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Open sidebar
OPAL
src
Commits
66f04040
Commit
66f04040
authored
Dec 10, 2020
by
ext-calvo_p
Browse files
Review cyclotron field map type
parent
3536b22f
Changes
6
Expand all
Hide whitespace changes
Inline
Side-by-side
Showing
6 changed files
with
226 additions
and
287 deletions
+226
-287
src/Algorithms/ParallelCyclotronTracker.cpp
src/Algorithms/ParallelCyclotronTracker.cpp
+16
-33
src/Classic/AbsBeamline/Cyclotron.cpp
src/Classic/AbsBeamline/Cyclotron.cpp
+159
-200
src/Classic/AbsBeamline/Cyclotron.h
src/Classic/AbsBeamline/Cyclotron.h
+41
-44
src/Distribution/ClosedOrbitFinder.h
src/Distribution/ClosedOrbitFinder.h
+2
-3
src/Distribution/Distribution.cpp
src/Distribution/Distribution.cpp
+4
-3
src/Elements/OpalCyclotron.cpp
src/Elements/OpalCyclotron.cpp
+4
-4
No files found.
src/Algorithms/ParallelCyclotronTracker.cpp
View file @
66f04040
...
...
@@ -396,8 +396,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
// useful information
spiral_flag
=
cycl_m
->
getSpiralFlag
();
if
(
spiral_flag
)
{
if
(
spiral_flag
)
{
*
gmsg
<<
endl
<<
"* This is a Spiral Inflector Simulation! This means the following:"
<<
endl
;
*
gmsg
<<
"* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!"
<<
endl
;
*
gmsg
<<
"* (Use BANDRF type cyclotron and use RFMAPFN to load both magnetic"
<<
endl
;
...
...
@@ -416,7 +415,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
}
// Fresh run (no restart):
if
(
!
OpalData
::
getInstance
()
->
inRestartRun
())
{
if
(
!
OpalData
::
getInstance
()
->
inRestartRun
())
{
// Get reference values from cyclotron element
// For now, these are still stored in mm. should be the only ones. -DW
...
...
@@ -426,7 +425,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
referencePr
=
cycl_m
->
getPRinit
();
referencePz
=
cycl_m
->
getPZinit
();
if
(
referenceTheta
<=
-
180.0
||
referenceTheta
>
180.0
)
{
if
(
referenceTheta
<=
-
180.0
||
referenceTheta
>
180.0
)
{
throw
OpalException
(
"Error in ParallelCyclotronTracker::visitCyclotron"
,
"PHIINIT is out of [-180, 180)!"
);
}
...
...
@@ -437,25 +436,22 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
float
insqrt
=
referencePtot
*
referencePtot
-
\
referencePr
*
referencePr
-
referencePz
*
referencePz
;
if
(
insqrt
<
0
)
{
if
(
insqrt
>
-
1.0e-10
)
{
if
(
insqrt
<
0
)
{
if
(
insqrt
>
-
1.0e-10
)
{
referencePt
=
0.0
;
}
else
{
throw
OpalException
(
"Error in ParallelCyclotronTracker::visitCyclotron"
,
"Pt imaginary!"
);
}
}
else
{
referencePt
=
std
::
sqrt
(
insqrt
);
}
if
(
referencePtot
<
0.0
)
if
(
referencePtot
<
0.0
)
{
referencePt
*=
-
1.0
;
}
// End calculate referencePt
// Restart a run:
...
...
@@ -463,15 +459,15 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
// If the user wants to save the restarted run in local frame,
// make sure the previous h5 file was local too
if
(
Options
::
psDumpFrame
!=
Options
::
GLOBAL
)
{
if
(
!
previousH5Local
)
{
if
(
Options
::
psDumpFrame
!=
Options
::
GLOBAL
)
{
if
(
!
previousH5Local
)
{
throw
OpalException
(
"Error in ParallelCyclotronTracker::visitCyclotron"
,
"You are trying a local restart from a global h5 file!"
);
}
}
// Else, if the user wants to save the restarted run in global frame,
// make sure the previous h5 file was global too
}
else
{
if
(
previousH5Local
)
{
}
else
{
if
(
previousH5Local
)
{
throw
OpalException
(
"Error in ParallelCyclotronTracker::visitCyclotron"
,
"You are trying a global restart from a local h5 file!"
);
}
...
...
@@ -481,7 +477,7 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
referencePhi
*=
Physics
::
deg2rad
;
referencePsi
*=
Physics
::
deg2rad
;
referencePtot
=
bega
;
if
(
referenceTheta
<=
-
180.0
||
referenceTheta
>
180.0
)
{
if
(
referenceTheta
<=
-
180.0
||
referenceTheta
>
180.0
)
{
throw
OpalException
(
"Error in ParallelCyclotronTracker::visitCyclotron"
,
"PHIINIT is out of [-180, 180)!"
);
}
...
...
@@ -539,21 +535,8 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
*
gmsg
<<
std
::
boolalpha
<<
"* Superpose electric field maps -> "
<<
superpose
<<
endl
;
}
/**
* To ease the initialise() function, set a integral parameter fieldflag internally.
* Its value is by the option "TYPE" of the element "CYCLOTRON"
* fieldflag = 1, readin PSI format measured field file (default)
* fieldflag = 2, readin carbon cyclotron field file created by Jianjun Yang, TYPE=CARBONCYCL
* fieldflag = 3, readin ANSYS format file for CYCIAE-100 created by Jianjun Yang, TYPE=CYCIAE
* fieldflag = 4, readin AVFEQ format file for Riken cyclotrons
* fieldflag = 5, readin FFA format file for MSU/FNAL FFA
* fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
* fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
*/
int
fieldflag
=
cycl_m
->
getFieldFlag
(
type
);
// Read in cyclotron field maps (midplane + 3D fields if desired).
cycl_m
->
initialise
(
itsBunch_m
,
fieldflag
,
cycl_m
->
getBScale
());
// Read in cyclotron field maps
cycl_m
->
initialise
(
itsBunch_m
,
cycl_m
->
getBScale
());
double
BcParameter
[
8
]
=
{};
...
...
@@ -3520,4 +3503,4 @@ void ParallelCyclotronTracker::initPathLength() {
// we need to reset the path length of each bunch
itsDataSink
->
setMultiBunchInitialPathLengh
(
mbHandler_m
.
get
());
}
}
\ No newline at end of file
}
src/Classic/AbsBeamline/Cyclotron.cpp
View file @
66f04040
This diff is collapsed.
Click to expand it.
src/Classic/AbsBeamline/Cyclotron.h
View file @
66f04040
...
...
@@ -35,7 +35,6 @@ class Fieldmap;
class
LossDataSink
;
class
TrimCoil
;
enum
BFieldType
{
PSIBF
,
CARBONBF
,
ANSYSBF
,
AVFEQBF
,
FFABF
,
BANDRF
,
SYNCHRO
};
struct
BfieldData
{
std
::
string
filename
;
...
...
@@ -60,19 +59,12 @@ struct BfieldData {
std
::
vector
<
double
>
g3
;
// for Btheta
// Grid-Size
//need to be read from inputfile.
int
nrad
,
ntet
;
// one more grid line is stored in azimuthal direction:
int
ntetS
;
// total grid points number.
int
ntot
;
int
nrad
,
ntet
;
//need to be read from inputfile.
int
ntetS
;
// one more grid line is stored in azimuthal direction
int
ntot
;
// total grid points number.
// Mean and Maximas
double
bacc
,
dbtmx
,
dbttmx
,
dbtttmx
;
};
struct
BPositions
{
...
...
@@ -88,16 +80,22 @@ struct BPositions {
};
// Class Cyclotron
// ------------------------------------------------------------------------
/// Interface for a Cyclotron.
// This class defines the abstract interface for a Cyclotron.
class
Cyclotron
:
public
Component
{
public:
enum
class
BFieldType
{
PSIBF
,
CARBONBF
,
ANSYSBF
,
AVFEQBF
,
FFABF
,
BANDRF
,
SYNCHRO
};
/// Constructor with given name.
explicit
Cyclotron
(
const
std
::
string
&
name
);
explicit
Cyclotron
(
const
std
::
string
&
name
);
Cyclotron
();
Cyclotron
(
const
Cyclotron
&
);
...
...
@@ -122,13 +120,11 @@ public:
void
setRFFCoeffFN
(
std
::
vector
<
std
::
string
>
rff_coeff_fn
);
void
setRFVCoeffFN
(
std
::
vector
<
std
::
string
>
rfv_coeff_fn
);
int
getFieldFlag
(
const
std
::
string
&
type
)
const
;
void
setType
(
std
::
string
t
);
void
setCyclotronType
(
std
::
string
t
);
const
std
::
string
&
getCyclotronType
()
const
;
virtual
ElementBase
::
ElementType
getType
()
const
;
virtual
void
getDimensions
(
double
&
zBegin
,
double
&
zEnd
)
const
;
virtual
void
getDimensions
(
double
&
zBegin
,
double
&
zEnd
)
const
;
unsigned
int
getNumberOfTrimcoils
()
const
;
...
...
@@ -165,7 +161,7 @@ public:
void
setEScale
(
std
::
vector
<
double
>
bs
);
virtual
double
getEScale
(
unsigned
int
i
)
const
;
void
setTrimCoils
(
const
std
::
vector
<
TrimCoil
*>
&
trimcoils
);
void
setTrimCoils
(
const
std
::
vector
<
TrimCoil
*>&
trimcoils
);
void
setSuperpose
(
std
::
vector
<
bool
>
flag
);
virtual
bool
getSuperpose
(
unsigned
int
i
)
const
;
...
...
@@ -191,17 +187,17 @@ public:
void
setSpiralFlag
(
bool
spiral_flag
);
virtual
bool
getSpiralFlag
()
const
;
virtual
bool
apply
(
const
size_t
&
id
,
const
double
&
t
,
Vector_t
&
E
,
Vector_t
&
B
);
virtual
bool
apply
(
const
size_t
&
id
,
const
double
&
t
,
Vector_t
&
E
,
Vector_t
&
B
);
virtual
bool
apply
(
const
Vector_t
&
R
,
const
Vector_t
&
P
,
const
double
&
t
,
Vector_t
&
E
,
Vector_t
&
B
);
virtual
bool
apply
(
const
Vector_t
&
R
,
const
Vector_t
&
P
,
const
double
&
t
,
Vector_t
&
E
,
Vector_t
&
B
);
virtual
void
apply
(
const
double
&
rad
,
const
double
&
z
,
const
double
&
tet_rad
,
double
&
br
,
double
&
bt
,
double
&
bz
);
virtual
void
initialise
(
PartBunchBase
<
double
,
3
>
*
bunch
,
double
&
startField
,
double
&
endField
);
virtual
void
initialise
(
PartBunchBase
<
double
,
3
>*
bunch
,
double
&
startField
,
double
&
endField
);
virtual
void
initialise
(
PartBunchBase
<
double
,
3
>
*
bunch
,
const
int
&
fieldflag
,
const
double
&
scaleFactor
);
virtual
void
initialise
(
PartBunchBase
<
double
,
3
>*
bunch
,
const
double
&
scaleFactor
);
virtual
void
finalise
();
...
...
@@ -217,34 +213,38 @@ public:
double
&
bt
,
double
&
bz
);
void
read
(
const
int
&
fieldflag
,
const
double
&
scaleFactor
);
void
read
(
const
double
&
scaleFactor
);
void
setBFieldType
();
BFieldType
fieldType_m
;
private:
/// Apply trim coils (calculate field contributions) with smooth field transition
void
applyTrimCoil
(
const
double
r
,
const
double
z
,
const
double
tet_rad
,
double
&
br
,
double
&
bz
);
/// Apply trim coils (calculate field contributions)
void
applyTrimCoil_m
(
const
double
r
,
const
double
z
,
const
double
tet_rad
,
double
*
br
,
double
*
bz
);
void
applyTrimCoil_m
(
const
double
r
,
const
double
z
,
const
double
tet_rad
,
double
*
br
,
double
*
bz
);
protected:
void
getdiffs
();
double
gutdf5d
(
double
*
f
,
double
dx
,
const
int
kor
,
const
int
krl
,
const
int
lpr
);
double
gutdf5d
(
double
*
f
,
double
dx
,
const
int
kor
,
const
int
krl
,
const
int
lpr
);
void
initR
(
double
rmin
,
double
dr
,
int
nrad
);
void
getFieldFromFile
(
const
double
&
scaleFactor
);
void
getFieldFromFile_Carbon
(
const
double
&
scaleFactor
);
void
getFieldFromFile_CYCIAE
(
const
double
&
scaleFactor
);
void
getFieldFromFile_AVFEQ
(
const
double
&
scaleFactor
);
void
getFieldFromFile_FFA
(
const
double
&
scaleFactor
);
void
getFieldFromFile_BandRF
(
const
double
&
scaleFactor
);
void
getFieldFromFile_Synchrocyclotron
(
const
double
&
scaleFactor
);
void
getFieldFromFile
_Ring
(
const
double
&
scaleFactor
);
void
getFieldFromFile_Carbon
(
const
double
&
scaleFactor
);
void
getFieldFromFile_CYCIAE
(
const
double
&
scaleFactor
);
void
getFieldFromFile_AVFEQ
(
const
double
&
scaleFactor
);
void
getFieldFromFile_FFA
(
const
double
&
scaleFactor
);
void
getFieldFromFile_BandRF
(
const
double
&
scaleFactor
);
void
getFieldFromFile_Synchrocyclotron
(
const
double
&
scaleFactor
);
inline
int
idx
(
int
irad
,
int
ktet
)
{
return
(
ktet
+
Bfield
.
ntetS
*
irad
);}
private:
std
::
string
fmapfn_m
;
/* stores the filename of the fieldmap */
...
...
@@ -267,7 +267,7 @@ private:
bool
spiral_flag_m
;
double
trimCoilThreshold_m
;
///< B-field threshold for applying trim coil
std
::
string
type_m
;
/
* what type of field we use */
std
::
string
type
Name
_m
;
/
/ name of the TYPE parameter in cyclotron
double
harm_m
;
double
bscale_m
;
// a scale factor for the B-field
...
...
@@ -287,9 +287,6 @@ private:
// Not implemented.
void
operator
=
(
const
Cyclotron
&
)
=
delete
;
BFieldType
myBFieldType_m
;
// RF field map handler
// Fieldmap *RFfield;
std
::
vector
<
Fieldmap
*>
RFfields_m
;
...
...
src/Distribution/ClosedOrbitFinder.h
View file @
66f04040
...
...
@@ -290,8 +290,7 @@ ClosedOrbitFinder<Value_type,
N_m
/=
cycl_m
->
getSymmetry
();
}
cycl_m
->
read
(
cycl_m
->
getFieldFlag
(
cycl_m
->
getCyclotronType
()),
cycl_m
->
getBScale
());
cycl_m
->
read
(
cycl_m
->
getBScale
());
// reserve storage for the orbit and momentum (--> size = 0, capacity = N_m+1)
/*
...
...
@@ -925,4 +924,4 @@ ClosedOrbitFinder<Value_type, Size_type, Stepper>::rotate(value_type angle, cons
}
#endif
\ No newline at end of file
#endif
src/Distribution/Distribution.cpp
View file @
66f04040
...
...
@@ -1225,9 +1225,10 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles,
else
*
gmsg
<<
"* SECTOR: "
<<
"match using single sector"
<<
endl
;
*
gmsg
<<
"* NSTEPS = "
<<
Nint
<<
endl
<<
"* HN= "
<<
CyclotronElement
->
getCyclHarm
()
<<
" PHIINIT= "
<<
CyclotronElement
->
getPHIinit
()
<<
endl
*
gmsg
<<
"* NSTEPS = "
<<
Nint
<<
endl
<<
"* HN = "
<<
CyclotronElement
->
getCyclHarm
()
<<
" PHIINIT = "
<<
CyclotronElement
->
getPHIinit
()
<<
endl
<<
"* FIELD MAP = "
<<
CyclotronElement
->
getFieldMapFN
()
<<
endl
<<
"* ----------------------------------------------------"
<<
endl
;
if
(
CyclotronElement
->
getFMLowE
()
<
0
||
...
...
src/Elements/OpalCyclotron.cpp
View file @
66f04040
...
...
@@ -170,7 +170,7 @@ void OpalCyclotron::update() {
cycl
->
setBScale
(
bscale
);
cycl
->
setType
(
type
);
cycl
->
set
Cyclotron
Type
(
type
);
cycl
->
setCyclHarm
(
harmnum
);
cycl
->
setMinR
(
minr
);
...
...
@@ -218,13 +218,13 @@ void OpalCyclotron::update() {
cycl
->
setRfFrequ
(
rff_str
);
cycl
->
setSuperpose
(
superpose_str
);
if
(
itsAttr
[
GEOMETRY
]
&&
obgeo_m
==
nullptr
)
{
if
(
itsAttr
[
GEOMETRY
]
&&
obgeo_m
==
nullptr
)
{
obgeo_m
=
BoundaryGeometry
::
find
(
Attributes
::
getString
(
itsAttr
[
GEOMETRY
]));
if
(
obgeo_m
)
{
if
(
obgeo_m
)
{
cycl
->
setBoundaryGeometry
(
obgeo_m
);
}
}
// Transmit "unknown" attributes.
OpalElement
::
updateUnknown
(
cycl
);
}
\ No newline at end of file
}
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment