diff --git a/CMakeLists.txt b/CMakeLists.txt
index cf7b3e703a563eb2bedfa077ce03b4b12fc3eae6..a8e1a29ee347ea15e58427de5117189bccb1aea5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -17,6 +17,15 @@ if (${CMAKE_BUILD_TYPE} STREQUAL "Release")
add_definitions (-DNOPAssert)
endif ()
+
+if (${CMAKE_BUILD_TYPE} STREQUAL "Debug")
+ add_compile_options (-fsanitize=undefined,address)
+ add_link_options (-fsanitize=undefined,address)
+endif ()
+
+
+
+
add_compile_options (-Wall)
add_compile_options (-Wunused)
add_compile_options (-Wextra)
@@ -56,7 +65,7 @@ message (STATUS "CMAKE_CXX_FLAGS_DEBUG: ${CMAKE_CXX_FLAGS_DEBUG}")
add_definitions (-DNOCTAssert)
#add_compile_options (-ferror-limit=1)
-#add_compile_options (-fsanitize=undefined,address)
+
#add_compile_options (-Wno-deprecated-declarations)
#add_compile_options (-Wno-unused)
#add_compile_options (-Wextra)
@@ -72,7 +81,7 @@ add_compile_options (-funroll-loops)
add_compile_options (-fstrict-aliasing)
add_compile_options (-DKOKKOS_DEPENDENCE)
add_compile_options (-Wno-return-type)
-add_compile_options (-gdwarf-2) # avoid dwarf errors on merlin
+add_compile_options (-gdwarf-4) # avoid dwarf errors on merlin
# Resolve all library dependencies
set (CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/CMakeModules")
diff --git a/input-files/test-0002.in b/input-files/test-0002.in
index de169a8c4c3150235d717c12a5587601ed3ebd2c..dee5559d8d2b65a800c36cf1ae2ef400ae9c36eb 100644
--- a/input-files/test-0002.in
+++ b/input-files/test-0002.in
@@ -17,7 +17,7 @@ Value,{OPALVERSION};
// Global Parameters
REAL rf_freq = 1.3e9; // RF frequency. (Hz)
-REAL n_particles = 1E6; // Number of particles in simulation.
+REAL n_particles = 1E5; // Number of particles in simulation.
REAL beam_bunch_charge = 1e-9; // Charge of bunch. (C)
// ----------------------------------------------------------------------------
@@ -32,6 +32,7 @@ REAL P0 = gamma*beta*EMASS; //inital z momentum
value, {Edes, P0};
D1: DRIFT, L = 1.0, ELEMEDGE = 0.0;
+myLine: Line = (D1);
BEAM1: BEAM, PARTICLE = ELECTRON, pc = P0, NPART = n_particles,
BFREQ = rf_freq, BCURRENT = beam_bunch_charge * rf_freq * 1E6, CHARGE = -1;
@@ -41,9 +42,9 @@ FS1: Fieldsolver, NX=32, NY=32, NZ=32, TYPE=FFT, PARFFTX = true, PARFFTY = true,
BCFFTX=OPEN, BCFFTY=OPEN, BCFFTZ=OPEN,
BBOXINCR = 1, GREENSF = INTEGRATED;
-Dist1: DISTRIBUTION, TYPE= GAUSS, SIGMAX=1, SIGMAY=1, SIGMAZ=1, SIGMAPX=1, SIGMAPY=1, SIGMAPZ=1;
+Dist1: DISTRIBUTION, TYPE= GAUSS, SIGMAX=1e-3, SIGMAY=1e-3, SIGMAZ=1e-3, SIGMAPX=1e-6, SIGMAPY=1e-6, SIGMAPZ=1e-6;
+
-myLine: Line = (D1);
TRACK, LINE = myLine, BEAM = BEAM1, MAXSTEPS = 1, DT = {1e-12}, ZSTOP=1.0;
RUN, METHOD = "PARALLEL", BEAM = BEAM1, FIELDSOLVER = FS1, DISTRIBUTION = Dist1;
diff --git a/input-files/test-0002.stat b/input-files/test-0002.stat
new file mode 100644
index 0000000000000000000000000000000000000000..6e876f49f9f18ba7171a89283b9f0b3f4aeb1252
--- /dev/null
+++ b/input-files/test-0002.stat
@@ -0,0 +1,304 @@
+SDDS1
+&description
+ text="Statistics data '../input-files/test-0002.in' 20/06/2024 23:23:03",
+ contents="stat parameters"
+&end
+¶meter
+ name=processors,
+ type=long,
+ description="Number of Cores used"
+&end
+¶meter
+ name=revision,
+ type=string,
+ description="git revision of opal"
+&end
+¶meter
+ name=flavor,
+ type=string,
+ description="OPAL flavor that wrote file"
+&end
+&column
+ name=t,
+ type=double,
+ units=ns,
+ description="1 Time"
+&end
+&column
+ name=s,
+ type=double,
+ units=m,
+ description="2 Path length"
+&end
+&column
+ name=numParticles,
+ type=long,
+ units=1,
+ description="3 Number of Macro Particles"
+&end
+&column
+ name=charge,
+ type=double,
+ units=1,
+ description="4 Bunch Charge"
+&end
+&column
+ name=energy,
+ type=double,
+ units=MeV,
+ description="5 Mean Bunch Energy"
+&end
+&column
+ name=rms_x,
+ type=double,
+ units=m,
+ description="6 RMS Beamsize in x"
+&end
+&column
+ name=rms_y,
+ type=double,
+ units=m,
+ description="7 RMS Beamsize in y"
+&end
+&column
+ name=rms_s,
+ type=double,
+ units=m,
+ description="8 RMS Beamsize in s"
+&end
+&column
+ name=rms_px,
+ type=double,
+ units=1,
+ description="9 RMS Normalized Momenta in x"
+&end
+&column
+ name=rms_py,
+ type=double,
+ units=1,
+ description="10 RMS Normalized Momenta in y"
+&end
+&column
+ name=rms_ps,
+ type=double,
+ units=1,
+ description="11 RMS Normalized Momenta in s"
+&end
+&column
+ name=emit_x,
+ type=double,
+ units=m,
+ description="12 Normalized Emittance x"
+&end
+&column
+ name=emit_y,
+ type=double,
+ units=m,
+ description="13 Normalized Emittance y"
+&end
+&column
+ name=emit_s,
+ type=double,
+ units=m,
+ description="14 Normalized Emittance s"
+&end
+&column
+ name=mean_x,
+ type=double,
+ units=m,
+ description="15 Mean Beam Position in x"
+&end
+&column
+ name=mean_y,
+ type=double,
+ units=m,
+ description="16 Mean Beam Position in y"
+&end
+&column
+ name=mean_s,
+ type=double,
+ units=m,
+ description="17 Mean Beam Position in s"
+&end
+&column
+ name=ref_x,
+ type=double,
+ units=m,
+ description="18 x coordinate of reference particle in lab cs"
+&end
+&column
+ name=ref_y,
+ type=double,
+ units=m,
+ description="19 y coordinate of reference particle in lab cs"
+&end
+&column
+ name=ref_z,
+ type=double,
+ units=m,
+ description="20 z coordinate of reference particle in lab cs"
+&end
+&column
+ name=ref_px,
+ type=double,
+ units=1,
+ description="21 x momentum of reference particle in lab cs"
+&end
+&column
+ name=ref_py,
+ type=double,
+ units=1,
+ description="22 y momentum of reference particle in lab cs"
+&end
+&column
+ name=ref_pz,
+ type=double,
+ units=1,
+ description="23 z momentum of reference particle in lab cs"
+&end
+&column
+ name=max_x,
+ type=double,
+ units=m,
+ description="24 Max Beamsize in x"
+&end
+&column
+ name=max_y,
+ type=double,
+ units=m,
+ description="25 Max Beamsize in y"
+&end
+&column
+ name=max_s,
+ type=double,
+ units=m,
+ description="26 Max Beamsize in s"
+&end
+&column
+ name=xpx,
+ type=double,
+ units=1,
+ description="27 Correlation xpx"
+&end
+&column
+ name=ypy,
+ type=double,
+ units=1,
+ description="28 Correlation ypy"
+&end
+&column
+ name=zpz,
+ type=double,
+ units=1,
+ description="29 Correlation zpz"
+&end
+&column
+ name=Dx,
+ type=double,
+ units=m,
+ description="30 Dispersion in x"
+&end
+&column
+ name=DDx,
+ type=double,
+ units=1,
+ description="31 Derivative of dispersion in x"
+&end
+&column
+ name=Dy,
+ type=double,
+ units=m,
+ description="32 Dispersion in y"
+&end
+&column
+ name=DDy,
+ type=double,
+ units=1,
+ description="33 Derivative of dispersion in y"
+&end
+&column
+ name=Bx_ref,
+ type=double,
+ units=T,
+ description="34 Bx-Field component of ref particle"
+&end
+&column
+ name=By_ref,
+ type=double,
+ units=T,
+ description="35 By-Field component of ref particle"
+&end
+&column
+ name=Bz_ref,
+ type=double,
+ units=T,
+ description="36 Bz-Field component of ref particle"
+&end
+&column
+ name=Ex_ref,
+ type=double,
+ units=MV/m,
+ description="37 Ex-Field component of ref particle"
+&end
+&column
+ name=Ey_ref,
+ type=double,
+ units=MV/m,
+ description="38 Ey-Field component of ref particle"
+&end
+&column
+ name=Ez_ref,
+ type=double,
+ units=MV/m,
+ description="39 Ez-Field component of ref particle"
+&end
+&column
+ name=dE,
+ type=double,
+ units=MeV,
+ description="40 energy spread of the beam"
+&end
+&column
+ name=dt,
+ type=double,
+ units=ns,
+ description="41 time step size"
+&end
+&column
+ name=partsOutside,
+ type=double,
+ units=1,
+ description="42 outside n*sigma of the beam"
+&end
+&column
+ name=DebyeLength,
+ type=double,
+ units=m,
+ description="43 Debye length in the boosted frame"
+&end
+&column
+ name=plasmaParameter,
+ type=double,
+ units=1,
+ description="44 Plasma parameter that gives no. of particles in a Debye sphere"
+&end
+&column
+ name=temperature,
+ type=double,
+ units=K,
+ description="45 Temperature of the beam"
+&end
+&column
+ name=rmsDensity,
+ type=double,
+ units=1,
+ description="46 RMS number density of the beam"
+&end
+&data
+ mode=ascii,
+ no_row_counts=1
+&end
+1
+OPALX 2024.1.00 git rev. #640dbb1f34e2dc23314bad26c94148a0e578fe17
+opal-t
+1.000000000000000e-03 2.997924188990765e-04 100000 -1.000000000000000e-09 9.999999999999980e+02 9.886771617229491e-04 9.900456978581636e-04 9.861933284523727e-04 9.970719439431197e-07 9.971772595976683e-07 1.003251673042217e-06 9.857814897437524e-10 9.872476063403890e-10 9.893927462964761e-10 1.695475357332832e-20 1.706100538623190e-20 4.746420270707041e-20 1.695475356915250e-20 1.706100538100345e-20 2.997924188990765e-04 -2.727234331922033e-23 -3.414707447737305e-23 1.957950928190183e+03 3.000478311466378e-03 2.999914659247034e-03 2.996972298154323e-03 -1.249744196182289e-03 2.643507480662082e-03 3.857264382145191e-03 -7.000706356679798e-13 4.605817787904672e-15 3.761377802323194e-12 -1.333863982837420e-16 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.621166691380730e-04 1.000000000000000e-03 0 3.883428884728593e-10 7.942613332657146e-15 1.025293560371121e-09 3.237641375683589e+13
diff --git a/input-files/test-0003.in b/input-files/test-0003.in
index 404458675b6a10710d8f343ef4618a84f2d13963..38289e38d25b45619bfb069c617d6a0b5b456f79 100644
--- a/input-files/test-0003.in
+++ b/input-files/test-0003.in
@@ -1,5 +1,5 @@
/*
- Test simple Gaussian distribution in a drift with space charge
+ Test simple Gaussian distribution
*/
OPTION, PSDUMPFREQ = 10; // 6d data written every 10th time step (h5).
@@ -9,7 +9,7 @@ OPTION, AUTOPHASE = 4; // Autophase is on, and phase of max energy
// gain will be found automatically for cavities.
OPTION, VERSION=10900;
-Title, string="Test simple Gaussian distribution ";
+Title, string="Test simple Gaussian distribution with space charge";
Value,{OPALVERSION};
@@ -32,6 +32,7 @@ REAL P0 = gamma*beta*EMASS; //inital z momentum
value, {Edes, P0};
D1: DRIFT, L = 1.0, ELEMEDGE = 0.0;
+myLine: Line = (D1);
BEAM1: BEAM, PARTICLE = ELECTRON, pc = P0, NPART = n_particles,
BFREQ = rf_freq, BCURRENT = beam_bunch_charge * rf_freq * 1E6, CHARGE = -1;
@@ -41,11 +42,11 @@ FS1: Fieldsolver, NX=32, NY=32, NZ=32, TYPE=FFT, PARFFTX = true, PARFFTY = true,
BCFFTX=OPEN, BCFFTY=OPEN, BCFFTZ=OPEN,
BBOXINCR = 1, GREENSF = INTEGRATED;
-Dist1: DISTRIBUTION, TYPE= GAUSS, SIGMAX=1, SIGMAY=1, SIGMAZ=1, SIGMAPX=1, SIGMAPY=1, SIGMAPZ=1;
+Dist1: DISTRIBUTION, TYPE= GAUSS, SIGMAX=1e-3, SIGMAY=1e-3, SIGMAZ=1e-3, SIGMAPX=1e-6, SIGMAPY=1e-6, SIGMAPZ=1e-6;
+
-myLine: Line = (D1);
-TRACK, LINE = myLine, BEAM = BEAM1, MAXSTEPS = 10, DT = {1e-12}, ZSTOP=1.0;
+TRACK, LINE = myLine, BEAM = BEAM1, MAXSTEPS = 2, DT = {1e-12}, ZSTOP=1.0;
RUN, METHOD = "PARALLEL", BEAM = BEAM1, FIELDSOLVER = FS1, DISTRIBUTION = Dist1;
ENDTRACK;
QUIT;
diff --git a/input-files/test-0003.stat b/input-files/test-0003.stat
new file mode 100644
index 0000000000000000000000000000000000000000..6cb8c4d2c9942835ae03716d5fbcb7e41f35d949
--- /dev/null
+++ b/input-files/test-0003.stat
@@ -0,0 +1,305 @@
+SDDS1
+&description
+ text="Statistics data '../input-files/test-0003.in' 20/06/2024 23:34:58",
+ contents="stat parameters"
+&end
+¶meter
+ name=processors,
+ type=long,
+ description="Number of Cores used"
+&end
+¶meter
+ name=revision,
+ type=string,
+ description="git revision of opal"
+&end
+¶meter
+ name=flavor,
+ type=string,
+ description="OPAL flavor that wrote file"
+&end
+&column
+ name=t,
+ type=double,
+ units=ns,
+ description="1 Time"
+&end
+&column
+ name=s,
+ type=double,
+ units=m,
+ description="2 Path length"
+&end
+&column
+ name=numParticles,
+ type=long,
+ units=1,
+ description="3 Number of Macro Particles"
+&end
+&column
+ name=charge,
+ type=double,
+ units=1,
+ description="4 Bunch Charge"
+&end
+&column
+ name=energy,
+ type=double,
+ units=MeV,
+ description="5 Mean Bunch Energy"
+&end
+&column
+ name=rms_x,
+ type=double,
+ units=m,
+ description="6 RMS Beamsize in x"
+&end
+&column
+ name=rms_y,
+ type=double,
+ units=m,
+ description="7 RMS Beamsize in y"
+&end
+&column
+ name=rms_s,
+ type=double,
+ units=m,
+ description="8 RMS Beamsize in s"
+&end
+&column
+ name=rms_px,
+ type=double,
+ units=1,
+ description="9 RMS Normalized Momenta in x"
+&end
+&column
+ name=rms_py,
+ type=double,
+ units=1,
+ description="10 RMS Normalized Momenta in y"
+&end
+&column
+ name=rms_ps,
+ type=double,
+ units=1,
+ description="11 RMS Normalized Momenta in s"
+&end
+&column
+ name=emit_x,
+ type=double,
+ units=m,
+ description="12 Normalized Emittance x"
+&end
+&column
+ name=emit_y,
+ type=double,
+ units=m,
+ description="13 Normalized Emittance y"
+&end
+&column
+ name=emit_s,
+ type=double,
+ units=m,
+ description="14 Normalized Emittance s"
+&end
+&column
+ name=mean_x,
+ type=double,
+ units=m,
+ description="15 Mean Beam Position in x"
+&end
+&column
+ name=mean_y,
+ type=double,
+ units=m,
+ description="16 Mean Beam Position in y"
+&end
+&column
+ name=mean_s,
+ type=double,
+ units=m,
+ description="17 Mean Beam Position in s"
+&end
+&column
+ name=ref_x,
+ type=double,
+ units=m,
+ description="18 x coordinate of reference particle in lab cs"
+&end
+&column
+ name=ref_y,
+ type=double,
+ units=m,
+ description="19 y coordinate of reference particle in lab cs"
+&end
+&column
+ name=ref_z,
+ type=double,
+ units=m,
+ description="20 z coordinate of reference particle in lab cs"
+&end
+&column
+ name=ref_px,
+ type=double,
+ units=1,
+ description="21 x momentum of reference particle in lab cs"
+&end
+&column
+ name=ref_py,
+ type=double,
+ units=1,
+ description="22 y momentum of reference particle in lab cs"
+&end
+&column
+ name=ref_pz,
+ type=double,
+ units=1,
+ description="23 z momentum of reference particle in lab cs"
+&end
+&column
+ name=max_x,
+ type=double,
+ units=m,
+ description="24 Max Beamsize in x"
+&end
+&column
+ name=max_y,
+ type=double,
+ units=m,
+ description="25 Max Beamsize in y"
+&end
+&column
+ name=max_s,
+ type=double,
+ units=m,
+ description="26 Max Beamsize in s"
+&end
+&column
+ name=xpx,
+ type=double,
+ units=1,
+ description="27 Correlation xpx"
+&end
+&column
+ name=ypy,
+ type=double,
+ units=1,
+ description="28 Correlation ypy"
+&end
+&column
+ name=zpz,
+ type=double,
+ units=1,
+ description="29 Correlation zpz"
+&end
+&column
+ name=Dx,
+ type=double,
+ units=m,
+ description="30 Dispersion in x"
+&end
+&column
+ name=DDx,
+ type=double,
+ units=1,
+ description="31 Derivative of dispersion in x"
+&end
+&column
+ name=Dy,
+ type=double,
+ units=m,
+ description="32 Dispersion in y"
+&end
+&column
+ name=DDy,
+ type=double,
+ units=1,
+ description="33 Derivative of dispersion in y"
+&end
+&column
+ name=Bx_ref,
+ type=double,
+ units=T,
+ description="34 Bx-Field component of ref particle"
+&end
+&column
+ name=By_ref,
+ type=double,
+ units=T,
+ description="35 By-Field component of ref particle"
+&end
+&column
+ name=Bz_ref,
+ type=double,
+ units=T,
+ description="36 Bz-Field component of ref particle"
+&end
+&column
+ name=Ex_ref,
+ type=double,
+ units=MV/m,
+ description="37 Ex-Field component of ref particle"
+&end
+&column
+ name=Ey_ref,
+ type=double,
+ units=MV/m,
+ description="38 Ey-Field component of ref particle"
+&end
+&column
+ name=Ez_ref,
+ type=double,
+ units=MV/m,
+ description="39 Ez-Field component of ref particle"
+&end
+&column
+ name=dE,
+ type=double,
+ units=MeV,
+ description="40 energy spread of the beam"
+&end
+&column
+ name=dt,
+ type=double,
+ units=ns,
+ description="41 time step size"
+&end
+&column
+ name=partsOutside,
+ type=double,
+ units=1,
+ description="42 outside n*sigma of the beam"
+&end
+&column
+ name=DebyeLength,
+ type=double,
+ units=m,
+ description="43 Debye length in the boosted frame"
+&end
+&column
+ name=plasmaParameter,
+ type=double,
+ units=1,
+ description="44 Plasma parameter that gives no. of particles in a Debye sphere"
+&end
+&column
+ name=temperature,
+ type=double,
+ units=K,
+ description="45 Temperature of the beam"
+&end
+&column
+ name=rmsDensity,
+ type=double,
+ units=1,
+ description="46 RMS number density of the beam"
+&end
+&data
+ mode=ascii,
+ no_row_counts=1
+&end
+1
+OPALX 2024.1.00 git rev. #640dbb1f34e2dc23314bad26c94148a0e578fe17
+opal-t
+1.000000000000000e-03 2.997924188990765e-04 100000 -1.000000000000000e-09 9.999999999999980e+02 9.886771617229491e-04 9.900456978581636e-04 9.861933284523727e-04 9.970719439431197e-07 9.971772595976683e-07 1.003251673042217e-06 9.857814897437524e-10 9.872476063403890e-10 9.893927462964761e-10 1.695475357332832e-20 1.706100538623190e-20 4.746420270707041e-20 1.695475356915250e-20 1.706100538100345e-20 2.997924188990765e-04 -2.727234331922033e-23 -3.414707447737305e-23 1.957950928190183e+03 3.000478311466378e-03 2.999914659247034e-03 2.996972298154323e-03 -1.249744196182289e-03 2.643507480662082e-03 3.857264382145191e-03 -7.000706356679798e-13 4.605817787904672e-15 3.761377802323194e-12 -1.333863982837420e-16 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.621166691380730e-04 1.000000000000000e-03 0 3.883428884728593e-10 7.942613332657146e-15 1.025293560371121e-09 3.237641375683589e+13
+2.000000000000000e-03 8.993772566972294e-04 100000 -1.000000000000000e-09 9.999999999999980e+02 9.886771617229491e-04 9.900456978581636e-04 9.861933284523727e-04 9.970719439431197e-07 9.971772595976683e-07 1.003251673042217e-06 9.857814897437524e-10 9.872476063403890e-10 9.893927462964761e-10 1.695475357332832e-20 1.706100538623190e-20 4.746420270707041e-20 1.695475356497669e-20 1.706100537577501e-20 5.995848377981529e-04 -2.727234331922033e-23 -3.414707447737305e-23 1.957950928190183e+03 3.000478311466378e-03 2.999914659247034e-03 2.996972298154323e-03 -1.249744196182289e-03 2.643507480662082e-03 3.857264382145191e-03 -7.000706356679798e-13 4.605817787904672e-15 3.761377802323194e-12 -1.333863982837420e-16 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.621166691380730e-04 1.000000000000000e-03 0 3.883407084178780e-10 7.942568744911947e-15 1.025293560371121e-09 3.237677726526380e+13
diff --git a/src/AbsBeamline/Probe.cpp b/src/AbsBeamline/Probe.cpp
index 0304b6cf9d087070acf13b3d1d9121708ec7dc6c..12dab2d04e22fec64184802a2479dc459bb3019b 100644
--- a/src/AbsBeamline/Probe.cpp
+++ b/src/AbsBeamline/Probe.cpp
@@ -99,8 +99,8 @@ bool Probe::doCheck(PartBunch_t* bunch, const int turnnumber, const double t, co
// dist1 > 0, right hand, dt > 0; dist1 < 0, left hand, dt < 0
double dist1 = (A_m * bunch->R(i)(0) + B_m * bunch->R(i)(1) + C_m) / R_m; // [m]
double dist2 = dist1 * std::sqrt(1.0 + 1.0 / tangle / tangle);
- double dt =
- dist2 / (std::sqrt(1.0 - 1.0 / (1.0 + dot(bunch->P(i), bunch->P(i)))) * Physics::c);
+ //double dt =
+ // dist2 / (std::sqrt(1.0 - 1.0 / (1.0 + dot(bunch->P(i), bunch->P(i)))) * Physics::c);
// ADA probepoint = bunch->R(i) + dist2 * bunch->P(i) / euclidean_norm(bunch->P(i));
probepoint = bunch->R(i) + dist2 * bunch->P(i) / std::sqrt(dot(bunch->P(i), bunch->P(i)));
diff --git a/src/AbsBeamline/Ring.cpp b/src/AbsBeamline/Ring.cpp
index 77e61b6a25542d618e8fe48b1ef5e142cc088e09..5c8c5d66426bc7674d4b4e0760f1b6a388b2f9a2 100644
--- a/src/AbsBeamline/Ring.cpp
+++ b/src/AbsBeamline/Ring.cpp
@@ -125,7 +125,7 @@ bool Ring::apply(
const double Q = Qview(id);
const double M = Mview(id);
- const short int Bin = Binview(id);
+ // const short int Bin = Binview(id);
gmsgALL << getName() << ": particle " << id << " at " << R
<< " m out of the field map boundary" << endl;
diff --git a/src/AbstractObjects/OpalData.h b/src/AbstractObjects/OpalData.h
index adffe2fe70058dae3009c78a907f16f44db78f9d..1260d0d8f7b32d647fabbac2d8eaf8d5ee748bb7 100644
--- a/src/AbstractObjects/OpalData.h
+++ b/src/AbstractObjects/OpalData.h
@@ -281,7 +281,7 @@ public:
private:
static bool isInstantiated;
- static OpalData* instance;
+ static OpalData* instance; // \todo should be a smart pointer and we then should get ridd of deleteInstance
static std::stack<OpalData*> stashedInstances;
OpalData();
diff --git a/src/Algorithms/ParallelTracker.cpp b/src/Algorithms/ParallelTracker.cpp
index 5601f5d5430d36fe91748e6e32a73bfed7a52517..e74b9d1088a5776e21f125b483e4b60cf3e08137 100644
--- a/src/Algorithms/ParallelTracker.cpp
+++ b/src/Algorithms/ParallelTracker.cpp
@@ -390,7 +390,6 @@ void ParallelTracker::execute() {
changeDT(back_track);
for (; step < trackSteps; ++step) {
- *gmsg << "at the beginning of loop " << endl;
Vector_t<double, 3> rmin(0.0), rmax(0.0);
if (itsBunch_m->getTotalNum() > 0) {
itsBunch_m->get_bounds(rmin, rmax);
@@ -442,7 +441,6 @@ void ParallelTracker::execute() {
if (std::abs(stepSizes_m.getZStop() - pathLength_m) < 0.5 * driftPerTimeStep) {
break;
}
- *gmsg << "at the end of loop pathLength_m= " << pathLength_m << endl;
}
if (globalEOL_m)
@@ -551,10 +549,12 @@ void ParallelTracker::emitParticles(long long step) {
}
void ParallelTracker::computeSpaceChargeFields(unsigned long long step) {
+
if (!itsBunch_m->hasFieldSolver()) {
+ *gmsg << "no solver avaidable " << endl;
return;
}
-
+
itsBunch_m->calcBeamParameters();
Quaternion alignment = getQuaternion(itsBunch_m->get_pmean(), Vector_t<double, 3>(0, 0, 1));
@@ -586,7 +586,9 @@ void ParallelTracker::computeSpaceChargeFields(unsigned long long step) {
h_Rot( i, j ) = rot(i, j);
}
}
+
Kokkos::deep_copy( Rot, h_Rot );
+
auto Rview = itsBunch_m->getParticleContainer()->R.getView();
auto Eview = itsBunch_m->getParticleContainer()->E.getView();
auto Bview = itsBunch_m->getParticleContainer()->B.getView();
diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp
index 568a870e32de16cb42f1a40b2355a01b57034b96..bceabbde8f2285289cc1253f9a4287f1b247cd3c 100644
--- a/src/Distribution/Distribution.cpp
+++ b/src/Distribution/Distribution.cpp
@@ -63,6 +63,18 @@ namespace DISTRIBUTION {
enum { TYPE, FNAME, SIGMAX, SIGMAY, SIGMAZ, SIGMAPX, SIGMAPY, SIGMAPZ, SIZE, CUTOFFPX, CUTOFFPY, CUTOFFPZ, CUTOFFX, CUTOFFY, CUTOFFLONG };
}
+/*
+namespace {
+ matrix_t getUnit6x6() {
+ matrix_t unit6x6(6, 6, 0.0); // Initialize a 6x6 matrix with all elements as 0.0
+ for (unsigned int i = 0; i < 6u; ++i) {
+ unit6x6(i, i) = 1.0; // Set diagonal elements to 1.0
+ }
+ return unit6x6;
+ }
+}
+*/
+
Distribution::Distribution()
: Definition(
DISTRIBUTION::SIZE, "DISTRIBUTION",
diff --git a/src/Main.cpp b/src/Main.cpp
index 6a6203518d36f0a6e1a865d44ae1f75bf0e5cacf..9784cc7795447ed7d0f63bd4f972fb333905db70 100644
--- a/src/Main.cpp
+++ b/src/Main.cpp
@@ -137,7 +137,7 @@ namespace OPALXMAIN {
int main(int argc, char* argv[]) {
ippl::initialize(argc, argv);
{
- gmsg = new Inform("OPAL");
+ gmsg = new Inform("OPAL-X");
namespace fs = boost::filesystem;
H5SetVerbosityLevel(1); // 65535);
@@ -199,7 +199,8 @@ int main(int argc, char* argv[]) {
if (is) {
*gmsg << "Reading startup file '" << startup << "'" << endl;
- parser.run(is);
+ parser.
+ run(is);
*gmsg << "Finished reading startup file." << endl;
}
FileStream::setEcho(Options::echo);
@@ -402,7 +403,9 @@ int main(int argc, char* argv[]) {
} catch (EarlyLeaveException& ex) {
// do nothing here
- } catch (OpalException& ex) {
+ }
+
+ catch (OpalException& ex) {
Inform errorMsg("Error", std::cerr, INFORM_ALL_NODES);
errorMsg << "\n*** User error detected by function \"" << ex.where() << "\"\n";
// stat->printWhere(errorMsg, true);
@@ -512,10 +515,14 @@ int main(int argc, char* argv[]) {
IpplTimings::print(
std::string("timing.dat"), OpalData::getInstance()->getProblemCharacteristicValues());
+
ippl::Comm->barrier();
Fieldmap::clearDictionary();
- OpalData::deleteInstance();
+
+ // \todo we should not need this OpalData::deleteInstance();
+
delete gmsg;
+ ippl::finalize();
return 0;
}
}
diff --git a/src/PartBunch/#[-Wunused-variable]# b/src/PartBunch/#[-Wunused-variable]#
new file mode 100644
index 0000000000000000000000000000000000000000..091276e8d88504da43085e1b5da7541ad81073bd
--- /dev/null
+++ b/src/PartBunch/#[-Wunused-variable]#
@@ -0,0 +1,2 @@
+ 499 | int tag = 101;
+
\ No newline at end of file
diff --git a/src/PartBunch/CMakeLists.txt b/src/PartBunch/CMakeLists.txt
index 9c9ec923fecb2316911022e4bb5e5113617f7322..6b3d5af74e7714d92957765575b3476781265581 100644
--- a/src/PartBunch/CMakeLists.txt
+++ b/src/PartBunch/CMakeLists.txt
@@ -1,5 +1,6 @@
set (_SRCS
PartBunch.cpp
+ FieldSolver.cpp
)
include_directories (
diff --git a/src/PartBunch/FieldSolver.cpp b/src/PartBunch/FieldSolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..98c16de7301a345ee3a0af48d4e062e23a72e9f9
--- /dev/null
+++ b/src/PartBunch/FieldSolver.cpp
@@ -0,0 +1,140 @@
+#include "PartBunch/FieldSolver.hpp"
+
+extern Inform* gmsg;
+
+template <>
+void FieldSolver<double,3>::initOpenSolver() {
+ // if constexpr (Dim == 3) {
+ ippl::ParameterList sp;
+ sp.add("output_type", OpenSolver_t<double, 3>::GRAD);
+ sp.add("use_heffte_defaults", false);
+ sp.add("use_pencils", true);
+ sp.add("use_reorder", false);
+ sp.add("use_gpu_aware", true);
+ sp.add("comm", ippl::p2p_pl);
+ sp.add("r2c_direction", 0);
+ sp.add("algorithm", OpenSolver_t<double, 3>::HOCKNEY);
+ initSolverWithParams<OpenSolver_t<double, 3>>(sp);
+ // } else {
+ //throw std::runtime_error("Unsupported dimensionality for OPEN solver");
+ //}
+}
+
+template <>
+void FieldSolver<double,3>::initSolver() {
+ Inform m;
+ if (this->getStype() == "FFT") {
+ initOpenSolver();
+ } else if (this->getStype() == "CG") {
+ initCGSolver();
+ } else if (this->getStype() == "P3M") {
+ initP3MSolver();
+ } else if (this->getStype() == "FFTOPEN") {
+ initFFTSolver();
+ } else if (this->getStype() == "NONE") {
+ initNullSolver();
+ }
+ else {
+ m << "No solver matches the argument: " << this->getStype() << endl;
+ }
+}
+
+template <>
+void FieldSolver<double,3>::setPotentialBCs() {
+ // CG requires explicit periodic boundary conditions while the periodic Poisson solver
+ // simply assumes them
+ if (this->getStype() == "CG") {
+ typedef ippl::BConds<Field<double, 3>, 3> bc_type;
+ bc_type allPeriodic;
+ for (unsigned int i = 0; i < 2 * 3; ++i) {
+ allPeriodic[i] = std::make_shared<ippl::PeriodicFace<Field<double, 3>>>(i);
+ }
+ phi_m->setFieldBC(allPeriodic);
+ }
+ }
+
+template<>
+void FieldSolver<double,3>::runSolver() {
+ constexpr int Dim = 3;
+
+ if (this->getStype() == "CG") {
+ CGSolver_t<double, 3>& solver = std::get<CGSolver_t<double, 3>>(this->getSolver());
+ solver.solve();
+
+ if (ippl::Comm->rank() == 0) {
+ std::stringstream fname;
+ fname << "data/CG_";
+ fname << ippl::Comm->size();
+ fname << ".csv";
+
+ Inform log(NULL, fname.str().c_str(), Inform::APPEND);
+ int iterations = solver.getIterationCount();
+ // Assume the dummy solve is the first call
+ if (iterations == 0) {
+ log << "residue,iterations" << endl;
+ }
+ // Don't print the dummy solve
+ if (iterations > 0) {
+ log << solver.getResidue() << "," << iterations << endl;
+ }
+ }
+ ippl::Comm->barrier();
+ } else if (this->getStype() == "FFT") {
+ if constexpr (Dim == 2 || Dim == 3) {
+ std::get<OpenSolver_t<double, 3>>(this->getSolver()).solve();
+ }
+ } else if (this->getStype() == "P3M") {
+ if constexpr (Dim == 3) {
+ std::get<P3MSolver_t<double, 3>>(this->getSolver()).solve();
+ }
+ } else if (this->getStype() == "FFTOPEN") {
+ if constexpr (Dim == 3) {
+ std::get<FFTSolver_t<double, 3>>(this->getSolver()).solve();
+ }
+ } else {
+ throw std::runtime_error("Unknown solver type");
+ }
+}
+
+/*
+template<>
+void FieldSolver<double,3>::initNullSolver() {
+ // if constexpr (Dim == 2 || Dim == 3) {
+ //ippl::ParameterList sp;
+ throw std::runtime_error("Not implemented Null solver");
+ // } else {
+ //throw std::runtime_error("Unsupported dimensionality for Null solver");
+ //}
+}
+*/
+
+/*
+template <>
+void FieldSolver<double,3>::initCGSolver() {
+ ippl::ParameterList sp;
+ sp.add("output_type", CGSolver_t<double, 3>::GRAD);
+ // Increase tolerance in the 1D case
+ sp.add("tolerance", 1e-10);
+
+ initSolverWithParams<CGSolver_t<double, 3>>(sp);
+}
+
+template<>
+void FieldSolver<double,3>::initP3MSolver() {
+ // if constexpr (Dim == 3) {
+ ippl::ParameterList sp;
+ sp.add("output_type", P3MSolver_t<double, 3>::GRAD);
+ sp.add("use_heffte_defaults", false);
+ sp.add("use_pencils", true);
+ sp.add("use_reorder", false);
+ sp.add("use_gpu_aware", true);
+ sp.add("comm", ippl::p2p_pl);
+ sp.add("r2c_direction", 0);
+
+ initSolverWithParams<P3MSolver_t<double, 3>>(sp);
+ // } else {
+ // throw std::runtime_error("Unsupported dimensionality for P3M solver");
+ // }
+}
+
+*/
\ No newline at end of file
diff --git a/src/PartBunch/FieldSolver.hpp b/src/PartBunch/FieldSolver.hpp
index 2ff8dfa4b3d123adc460568a686db569085ccb71..7573fc2f296a188707bd9ec90d9ebe62ee364dc9 100644
--- a/src/PartBunch/FieldSolver.hpp
+++ b/src/PartBunch/FieldSolver.hpp
@@ -70,77 +70,14 @@ public:
phi_m = phi;
}
- void initSolver() override {
- Inform m("solver ");
- if (this->getStype() == "FFT") {
- initFFTSolver();
- } else if (this->getStype() == "CG") {
- initCGSolver();
- } else if (this->getStype() == "P3M") {
- initP3MSolver();
- } else if (this->getStype() == "OPEN") {
- initOpenSolver();
- } else if (this->getStype() == "NONE") {
- initNullSolver();
- }
- else {
- m << "No solver matches the argument: " << this->getStype() << endl;
- }
- }
+ void initOpenSolver();
- void setPotentialBCs() {
- // CG requires explicit periodic boundary conditions while the periodic Poisson solver
- // simply assumes them
- if (this->getStype() == "CG") {
- typedef ippl::BConds<Field<T, Dim>, Dim> bc_type;
- bc_type allPeriodic;
- for (unsigned int i = 0; i < 2 * Dim; ++i) {
- allPeriodic[i] = std::make_shared<ippl::PeriodicFace<Field<T, Dim>>>(i);
- }
- phi_m->setFieldBC(allPeriodic);
- }
- }
+ void initSolver() override ;
- void runSolver() override {
- if (this->getStype() == "CG") {
- CGSolver_t<T, Dim>& solver = std::get<CGSolver_t<T, Dim>>(this->getSolver());
- solver.solve();
-
- if (ippl::Comm->rank() == 0) {
- std::stringstream fname;
- fname << "data/CG_";
- fname << ippl::Comm->size();
- fname << ".csv";
-
- Inform log(NULL, fname.str().c_str(), Inform::APPEND);
- int iterations = solver.getIterationCount();
- // Assume the dummy solve is the first call
- if (iterations == 0) {
- log << "residue,iterations" << endl;
- }
- // Don't print the dummy solve
- if (iterations > 0) {
- log << solver.getResidue() << "," << iterations << endl;
- }
- }
- ippl::Comm->barrier();
- } else if (this->getStype() == "FFT") {
- if constexpr (Dim == 2 || Dim == 3) {
- std::get<FFTSolver_t<T, Dim>>(this->getSolver()).solve();
- }
- } else if (this->getStype() == "P3M") {
- if constexpr (Dim == 3) {
- std::get<P3MSolver_t<T, Dim>>(this->getSolver()).solve();
- }
- } else if (this->getStype() == "OPEN") {
- if constexpr (Dim == 3) {
- std::get<OpenSolver_t<T, Dim>>(this->getSolver()).solve();
- }
- } else {
- throw std::runtime_error("Unknown solver type");
- }
- }
+ void setPotentialBCs();
+ void runSolver() override;
+
template <typename Solver>
void initSolverWithParams(const ippl::ParameterList& sp) {
this->getSolver().template emplace<Solver>();
@@ -162,74 +99,23 @@ public:
}
}
- void initNullSolver() {
- if constexpr (Dim == 2 || Dim == 3) {
- ippl::ParameterList sp;
- throw std::runtime_error("Not implemented Null solver");
- } else {
- throw std::runtime_error("Unsupported dimensionality for Null solver");
- }
- }
-
+ void initNullSolver() { }
+
void initFFTSolver() {
- if constexpr (Dim == 2 || Dim == 3) {
- ippl::ParameterList sp;
- sp.add("output_type", FFTSolver_t<T, Dim>::GRAD);
- sp.add("use_heffte_defaults", false);
- sp.add("use_pencils", true);
- sp.add("use_reorder", false);
- sp.add("use_gpu_aware", true);
- sp.add("comm", ippl::p2p_pl);
- sp.add("r2c_direction", 0);
-
- initSolverWithParams<FFTSolver_t<T, Dim>>(sp);
- } else {
- throw std::runtime_error("Unsupported dimensionality for FFT solver");
- }
+ ippl::ParameterList sp;
+ sp.add("output_type", FFTSolver_t<double, 3>::GRAD);
+ sp.add("use_heffte_defaults", false);
+ sp.add("use_pencils", true);
+ sp.add("use_reorder", false);
+ sp.add("use_gpu_aware", true);
+ sp.add("comm", ippl::p2p_pl);
+ sp.add("r2c_direction", 0);
+ initSolverWithParams<FFTSolver_t<double, 3>>(sp);
}
+
+ void initCGSolver() { }
- void initCGSolver() {
- ippl::ParameterList sp;
- sp.add("output_type", CGSolver_t<T, Dim>::GRAD);
- // Increase tolerance in the 1D case
- sp.add("tolerance", 1e-10);
-
- initSolverWithParams<CGSolver_t<T, Dim>>(sp);
- }
+ void initP3MSolver() { }
- void initP3MSolver() {
- if constexpr (Dim == 3) {
- ippl::ParameterList sp;
- sp.add("output_type", P3MSolver_t<T, Dim>::GRAD);
- sp.add("use_heffte_defaults", false);
- sp.add("use_pencils", true);
- sp.add("use_reorder", false);
- sp.add("use_gpu_aware", true);
- sp.add("comm", ippl::p2p_pl);
- sp.add("r2c_direction", 0);
-
- initSolverWithParams<P3MSolver_t<T, Dim>>(sp);
- } else {
- throw std::runtime_error("Unsupported dimensionality for P3M solver");
- }
- }
-
- void initOpenSolver() {
- if constexpr (Dim == 3) {
- ippl::ParameterList sp;
- sp.add("output_type", OpenSolver_t<T, Dim>::GRAD);
- sp.add("use_heffte_defaults", false);
- sp.add("use_pencils", true);
- sp.add("use_reorder", false);
- sp.add("use_gpu_aware", true);
- sp.add("comm", ippl::p2p_pl);
- sp.add("r2c_direction", 0);
- sp.add("algorithm", OpenSolver_t<T, Dim>::HOCKNEY);
-
- initSolverWithParams<OpenSolver_t<T, Dim>>(sp);
- } else {
- throw std::runtime_error("Unsupported dimensionality for OPEN solver");
- }
- }
};
#endif
diff --git a/src/PartBunch/PartBunch.cpp b/src/PartBunch/PartBunch.cpp
index a6e2132408d06eaa9da9487160ad02ebb21c616c..42d7dd89c51120d9e41aca05d2cbc8baa7328d4d 100644
--- a/src/PartBunch/PartBunch.cpp
+++ b/src/PartBunch/PartBunch.cpp
@@ -1,7 +1,100 @@
+ // some debug output ------------------------------------------------------------
+ //
+ /*
+ Kokkos::parallel_for("print q", ippl::getRangePolicy(Qview),
+ KOKKOS_LAMBDA(const int i) {
+ if (i<5){
+ double myQ = Qview(i);
+ std::cout << "qi= " << myQ << std::endl;
+ }
+ });
+ */
+
#include "PartBunch/PartBunch.hpp"
#include <boost/numeric/ublas/io.hpp>
#include "Utilities/Util.h"
+extern Inform* gmsg;
+
+template<>
+void PartBunch<double,3>::setSolver(std::string solver) {
+
+ if (this->solver_m != "")
+ *gmsg << "* Warning solver already initiated but overwrite ..." << endl;
+
+ this->solver_m = solver;
+
+ this->fcontainer_m->initializeFields(this->solver_m);
+
+ this->setFieldSolver(std::make_shared<FieldSolver_t>(
+ this->solver_m, &this->fcontainer_m->getRho(), &this->fcontainer_m->getE(),
+ &this->fcontainer_m->getPhi()));
+
+ this->fsolver_m->initSolver();
+
+ /// ADA we need to be able to set a load balancer when not having a field solver
+ this->setLoadBalancer(std::make_shared<LoadBalancer_t>(
+ this->lbt_m, this->fcontainer_m, this->pcontainer_m, this->fsolver_m));
+
+ *gmsg << "* Solver and Load Balancer set" << endl;
+}
+
+template<>
+void PartBunch<double,3>::spaceChargeEFieldCheck() {
+
+ Inform msg("EParticleStats");
+
+ auto pE_view = this->pcontainer_m->E.getView();
+
+ double avgE = 0.0;
+ double minEComponent = std::numeric_limits<double>::max();
+ double maxEComponent = std::numeric_limits<double>::min();
+ double minE = std::numeric_limits<double>::max();
+ double maxE = std::numeric_limits<double>::min();
+ double cc = getCouplingConstant();
+
+ int myRank = ippl::Comm->rank();
+ Kokkos::parallel_reduce(
+ "check e-field", this->getLocalNum(),
+ KOKKOS_LAMBDA(const int i, double& loc_avgE, double& loc_minEComponent,
+ double& loc_maxEComponent, double& loc_minE, double& loc_maxE) {
+ double EX = pE_view[i][0]*cc;
+ double EY = pE_view[i][1]*cc;
+ double EZ = pE_view[i][2]*cc;
+
+ double ENorm = Kokkos::sqrt(EX*EX + EY*EY + EZ*EZ);
+
+ loc_avgE += ENorm;
+
+ loc_minEComponent = EX < loc_minEComponent ? EX : loc_minEComponent;
+ loc_minEComponent = EY < loc_minEComponent ? EY : loc_minEComponent;
+ loc_minEComponent = EZ < loc_minEComponent ? EZ : loc_minEComponent;
+
+ loc_maxEComponent = EX > loc_maxEComponent ? EX : loc_maxEComponent;
+ loc_maxEComponent = EY > loc_maxEComponent ? EY : loc_maxEComponent;
+ loc_maxEComponent = EZ > loc_maxEComponent ? EZ : loc_maxEComponent;
+
+ loc_minE = ENorm < loc_minE ? ENorm : loc_minE;
+ loc_maxE = ENorm > loc_maxE ? ENorm : loc_maxE;
+ },
+ Kokkos::Sum<double>(avgE), Kokkos::Min<double>(minEComponent),
+ Kokkos::Max<double>(maxEComponent), Kokkos::Min<double>(minE),
+ Kokkos::Max<double>(maxE));
+
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &avgE, &avgE, 1, MPI_DOUBLE, MPI_SUM, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &minEComponent, &minEComponent, 1, MPI_DOUBLE, MPI_MIN, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &maxEComponent, &maxEComponent, 1, MPI_DOUBLE, MPI_MAX, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &minE, &minE, 1, MPI_DOUBLE, MPI_MIN, 0, ippl::Comm->getCommunicator());
+ MPI_Reduce(myRank == 0 ? MPI_IN_PLACE : &maxE, &maxE, 1, MPI_DOUBLE, MPI_MAX, 0, ippl::Comm->getCommunicator());
+
+ avgE /= this->getTotalNum();
+ msg << "avgENorm = " << avgE << endl;
+ msg << "minEComponent = " << minEComponent << endl;
+ msg << "maxEComponent = " << maxEComponent << endl;
+ msg << "minE = " << minE << endl;
+ msg << "maxE = " << maxE << endl;
+}
+
template<>
void PartBunch<double,3>::calcBeamParameters() {
// Usefull constants
@@ -97,12 +190,16 @@ void PartBunch<double,3>::calcBeamParameters() {
rmax_m(i) = rmax[i];
rmin_m(i) = rmin[i];
}
-
ippl::Comm->barrier();
-
-
}
+template<>
+void PartBunch<double,3>::pre_run() {
+ this->fcontainer_m->getRho() = 0.0;
+ this->fsolver_m->runSolver();
+ *gmsg << "* Solver warmup done" << endl;
+}
+
template<>
Inform& PartBunch<double,3>::print(Inform& os) {
// if (this->getLocalNum() != 0) { // to suppress Nans
@@ -149,18 +246,18 @@ Inform& PartBunch<double,3>::print(Inform& os) {
template <>
void PartBunch<double,3>::bunchUpdate() {
- /* \frief
+
+ /* \brief
1. calculates and set hr
2. do repartitioning
*/
auto *mesh = &this->fcontainer_m->getMesh();
- auto *FL = &this->fcontainer_m->getFL();
+ auto *FL = &this->fcontainer_m->getFL();
std::shared_ptr<ParticleContainer_t> pc = this->getParticleContainer();
pc->computeMinMaxR();
-
/*
This needs to go
auto Rview = this->getParticleContainer()->R.getView();
@@ -172,20 +269,21 @@ void PartBunch<double,3>::bunchUpdate() {
*/
/// \brief assume o < 0.0?
+
ippl::Vector<double, 3> o = pc->getMinR();
ippl::Vector<double, 3> e = pc->getMaxR();
ippl::Vector<double, 3> l = e - o;
hr_m = (1.0+this->OPALFieldSolver_m->getBoxIncr()/100.)*(l / this->nr_m);
mesh->setMeshSpacing(hr_m);
mesh->setOrigin(o-0.5*hr_m*this->OPALFieldSolver_m->getBoxIncr()/100.);
-
+
pc->getLayout().updateLayout(*FL, *mesh);
pc->update();
this->getFieldContainer()->setRMin(o);
this->getFieldContainer()->setRMax(e);
this->getFieldContainer()->setHr(hr_m);
-
+
/* old Tracker
this->calcBeamParameters();
@@ -202,65 +300,115 @@ void PartBunch<double,3>::bunchUpdate() {
this->getParticleContainer()->getLayout().updateLayout(*FL, *mesh);
this->getParticleContainer()->update();
*/
-/*
+
this->isFirstRepartition_m = true;
this->loadbalancer_m->initializeORB(FL, mesh);
- this->loadbalancer_m->repartition(FL, mesh, this->isFirstRepartition_m);
-*/
- this->updateMoments();
+ // \fixme with the OPEN solver repartion does not work.
+ // this->loadbalancer_m->repartition(FL, mesh, this->isFirstRepartition_m);
+ this->updateMoments();
}
+/*
+
+ auto rhoView = rho.getView();
+ using index_array_type_3d = typename ippl::RangePolicy<Dim>::index_array_type;
+ const index_array_type_3d& args3d = {8,8,8};
+ auto res3d = ippl::apply(rhoView,args3d);
+ m << "Type of variable: " << typeid(res3d).name() << " rho[8,8,8]= " << res3d << endl;
+
+ double qsum = 0.0;
+ for (int i=0; i<nr_m[0]; i++) {
+ for (int j=0; j<nr_m[1]; j++) {
+ for (int k=0; k<nr_m[2]; k++) {
+ const index_array_type_3d& args3d = {i,j,k};
+ qsum += ippl::apply(rhoView,args3d);
+ }
+ }
+ }
+
+ m << "sum(rho_m)= " << qsum << endl;
+
+
+ auto Qview = this->pcontainer_m->Q.getView();
+ using index_array_type_1d = typename ippl::RangePolicy<1>::index_array_type;
+ const index_array_type_1d& args1d = {0};
+ auto res1d = ippl::apply(Qview,args1d);
+ m << "Type of variable: " << typeid(res1d).name() << " value " << res1d << endl;
+
+
+
+ */
+
+
+
+// ADA
template <>
void PartBunch<double,3>::computeSelfFields() {
- static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("SolveTimer");
- IpplTimings::startTimer(SolveTimer);
-
- Field_t<3>& rho = this->fcontainer_m->getRho();
- ippl::ParticleAttrib<double>& Q = this->pcontainer_m->Q;
- typename Base::particle_position_type& R = this->pcontainer_m->R;
-
- rho = 0.0;
- Q = Q*this->pcontainer_m->dt;
- this->qi_m = this->qi_m * getdT();
- scatter(Q, rho, R);
- Q = Q/this->pcontainer_m->dt;
- this->qi_m = this->qi_m / getdT();
- rho = rho/getdT();
-
- double gammaz = this->pcontainer_m->getMeanGammaZ();
- double scaleFactor = 1;
- // double scaleFactor = Physics::c * getdT();
- //and get meshspacings in real units [m]
- Vector_t<double, 3> hr_scaled = hr_m * scaleFactor;
- hr_scaled[2] *= gammaz;
-
- double localDensity2 = 0, Density2=0;
- double tmp2 = 1 / hr_scaled[0] * 1 / hr_scaled[1] * 1 / hr_scaled[2];
-
- using index_array_type = typename ippl::RangePolicy<Dim>::index_array_type;
-
- //divide charge by a 'grid-cube' volume to get [C/m^3]
- rho = rho *tmp2;
-
- double Npoints = nr_m[0] * nr_m[1] * nr_m[2];
-
- auto rhoView = rho.getView();
- localDensity2 = 0.;
- ippl::parallel_reduce(
- "Density stats", ippl::getRangePolicy(rhoView),
- KOKKOS_LAMBDA(const index_array_type& args, double& den) {
- double val = ippl::apply(rhoView, args);
- den += Kokkos::pow(val, 2);
- },
- Kokkos::Sum<double>(localDensity2) );
- ippl::Comm->reduce(localDensity2, Density2, 1, std::plus<double>());
-
- rmsDensity_m = std::sqrt( (1.0 /Npoints) * Density2 / Physics::q_e / Physics::q_e );
-
- this->calcDebyeLength();
-
- IpplTimings::stopTimer(SolveTimer);
+ Inform m("computeSelfFields ");
+ static IpplTimings::TimerRef SolveTimer = IpplTimings::getTimer("SolveTimer");
+ IpplTimings::startTimer(SolveTimer);
+
+ Field_t<3>& rho = this->fcontainer_m->getRho();
+ auto rhoView = rho.getView();
+ ippl::ParticleAttrib<double>& Q = this->pcontainer_m->Q;
+ typename Base::particle_position_type& R = this->pcontainer_m->R;
+
+ rho = 0.0;
+
+ Q = Q*this->pcontainer_m->dt;
+ this->qi_m = this->qi_m * getdT();
+
+ scatter(Q, rho, R);
+
+ Q = Q/this->pcontainer_m->dt;
+ this->qi_m = this->qi_m / getdT();
+ rho = rho/getdT();
+
+ double gammaz = this->pcontainer_m->getMeanGammaZ();
+ double scaleFactor = 1;
+ // double scaleFactor = Physics::c * getdT();
+ //and get meshspacings in real units [m]
+
+ Vector_t<double, 3> hr_scaled = hr_m * scaleFactor;
+ hr_scaled[2] *= gammaz;
+
+
+ double tmp2 = 1 / hr_scaled[0] * 1 / hr_scaled[1] * 1 / hr_scaled[2];
+
+ using index_array_type = typename ippl::RangePolicy<Dim>::index_array_type;
+
+ //divide charge by a 'grid-cube' volume to get [C/m^3]
+ rho = rho *tmp2;
+
+ double Npoints = nr_m[0] * nr_m[1] * nr_m[2];
+
+ double localDensity2 = 0.0, Density2=0.0;
+ ippl::parallel_reduce(
+ "Density stats", ippl::getRangePolicy(rhoView),
+ KOKKOS_LAMBDA(const index_array_type& args, double& den) {
+ double val = ippl::apply(rhoView, args);
+ den += Kokkos::pow(val, 2);
+ },
+ Kokkos::Sum<double>(localDensity2) );
+
+ ippl::Comm->reduce(localDensity2, Density2, 1, std::plus<double>());
+
+ rmsDensity_m = std::sqrt( (1.0 /Npoints) * Density2 / Physics::q_e / Physics::q_e );
+
+ this->calcDebyeLength();
+
+ m << "gammaz = " << gammaz << endl;
+ m << "hr_scaled = " << hr_scaled << endl;
+ m << "coupling constant= " << getCouplingConstant() << endl;
+
+ this->fsolver_m->runSolver();
+
+ gather(this->pcontainer_m->E, this->fcontainer_m->getE(), this->pcontainer_m->R);
+
+ spaceChargeEFieldCheck();
+
+ IpplTimings::stopTimer(SolveTimer);
}
template <>
@@ -279,10 +427,14 @@ void PartBunch<double,3>::scatterCIC() {
scatter(*q, *rho, *R);
+ m << "gammz= " << this->pcontainer_m->getMeanP()[2] << endl;
+
double relError = std::fabs((Q - (*rho).sum()) / Q);
size_type TotalParticles = 0;
size_type localParticles = this->pcontainer_m->getLocalNum();
+ m << "computeSelfFields sum rho " << (*rho).sum() << endl;
+
ippl::Comm->reduce(localParticles, TotalParticles, 1, std::plus<size_type>());
if (ippl::Comm->rank() == 0) {
diff --git a/src/PartBunch/PartBunch.hpp b/src/PartBunch/PartBunch.hpp
index 4d3fb66d3293a82bdf12135fba65cd4194a5b55c..b4aade3e1077763b2ca485219417ce9cca01e3e5 100644
--- a/src/PartBunch/PartBunch.hpp
+++ b/src/PartBunch/PartBunch.hpp
@@ -62,6 +62,14 @@ diagnostics(): calculate statistics and maybe write tp h5 and stat files
#include "Algorithms/PartData.h"
+
+template <typename T>
+KOKKOS_INLINE_FUNCTION typename T::value_type L2Norm(T& x) {
+ return sqrt(dot(x, x).apply());
+}
+
+
+
using view_type = typename ippl::detail::ViewType<ippl::Vector<double, 3>, 1>::view_type;
@@ -195,6 +203,12 @@ private:
double couplingConstant_m;
+ /// step in a TRACK command
+ long long localTrackStep_m;
+
+ /// if multiple TRACK commands
+ long long globalTrackStep_m;
+
// unit state of PartBunch
// UnitState_t unit_state_m;
// UnitState_t stateOfLastBoundP_m;
@@ -286,11 +300,6 @@ public:
IpplTimings::startTimer(prerun);
pre_run();
IpplTimings::stopTimer(prerun);
-
- /*
- end wrmup
- */
-
}
void bunchUpdate();
@@ -304,39 +313,10 @@ public:
return this->pcontainer_m;
}
- void setSolver(std::string solver) {
- Inform m("setSolver ");
-
- if (this->solver_m != "")
- m << "Warning solver already initiated but overwrite ..." << endl;
-
- this->solver_m = solver;
-
- this->fcontainer_m->initializeFields(this->solver_m);
-
- this->setFieldSolver(std::make_shared<FieldSolver_t>(
- this->solver_m, &this->fcontainer_m->getRho(), &this->fcontainer_m->getE(),
- &this->fcontainer_m->getPhi()));
-
- this->fsolver_m->initSolver();
-
- /// ADA we need to be able to set a load balancer when not having a field solver
- this->setLoadBalancer(std::make_shared<LoadBalancer_t>(
- this->lbt_m, this->fcontainer_m, this->pcontainer_m, this->fsolver_m));
- }
-
- void pre_run() override {
- static IpplTimings::TimerRef DummySolveTimer = IpplTimings::getTimer("solveWarmup");
- IpplTimings::startTimer(DummySolveTimer);
-
- this->fcontainer_m->getRho() = 0.0;
-
- this->fsolver_m->runSolver();
-
- IpplTimings::stopTimer(DummySolveTimer);
-
- }
+ void setSolver(std::string solver);
+ void pre_run() override ;
+
public:
void updateMoments(){
this->pcontainer_m->updateMoments();
@@ -399,10 +379,11 @@ public:
*/
double getCouplingConstant() const {
- return 1.0;
+ return couplingConstant_m;
}
-
+
void setCouplingConstant(double c) {
+ couplingConstant_m = c;
}
void calcBeamParameters();
@@ -584,14 +565,14 @@ public:
void setTEmission(double t) {
}
double getTEmission() {
- return 0.0;
+ return 0.0;
}
bool weHaveBins() {
- return false;
+ return false;
}
// void setPBins(PartBins* pbin) {}
size_t emitParticles(double eZ) {
- return 0;
+ return 0;
}
void updateNumTotal() {
}
@@ -896,6 +877,10 @@ public:
double calculateAngle(double x, double y) {
return 0.0;
}
+
+ // Sanity check functions
+ void spaceChargeEFieldCheck();
+
};
/* \todo
diff --git a/src/PartBunch/ParticleContainer.hpp b/src/PartBunch/ParticleContainer.hpp
index 4c4ba8085b60debd2eed914af54ac4f8aecf21a1..b0b5169ca2ede0e2a49fe441b99071e4d39e66fe 100644
--- a/src/PartBunch/ParticleContainer.hpp
+++ b/src/PartBunch/ParticleContainer.hpp
@@ -175,6 +175,7 @@ public:
void computeDebyeLength(double density){
size_t Np = this->getTotalNum();
distMoments_m.computeDebyeLength(this->R.getView(), this->P.getView(), Np, density);
+ return distMoments_m.getDebyeLength();
}
private:
diff --git a/src/Structure/H5PartWrapper.cpp b/src/Structure/H5PartWrapper.cpp
index 212567b6133ca60e83d336c43c0a8d124406ae36..4180ea11f446caff2323b3a0ea479e78e5234cb9 100644
--- a/src/Structure/H5PartWrapper.cpp
+++ b/src/Structure/H5PartWrapper.cpp
@@ -484,7 +484,7 @@ void H5PartWrapper::copyStepData(h5_file_t source) {
void H5PartWrapper::sendFailureMessage(
bool failed, const std::string& where, const std::string& what) {
- int tag = 101;
+ // int tag = 101;
/* \todo Message* mess = new Message();
putMessage(*mess, failed);
Ippl::Comm->broadcast_all(mess, tag);
@@ -496,8 +496,8 @@ void H5PartWrapper::sendFailureMessage(
void H5PartWrapper::receiveFailureMessage(
int sourceNode, const std::string& where, const std::string& what) {
- int tag = 101;
- bool failed;
+ // int tag = 101;
+ bool failed=false;
/* \todo
Message* mess = Ippl::Comm->receive_block(sourceNode, tag);
getMessage(*mess, failed);
diff --git a/src/Structure/IpplInfoWrapper.cpp b/src/Structure/IpplInfoWrapper.cpp
index 2b236caec50f8cedde59969d1ea437afbddb2c98..6821b8da20dbbfde10f8e7dd06b34dcd62636d28 100644
--- a/src/Structure/IpplInfoWrapper.cpp
+++ b/src/Structure/IpplInfoWrapper.cpp
@@ -39,7 +39,7 @@ IpplInfoWrapper::IpplInfoWrapper(
arg_m[5] = buffer_m + warn_m;
arg_m[6] = buffer_m + warnLevel_m;
- int narg = 5;
+ // int narg = 5;
// instance_m = new Ippl(narg, arg_m, Ippl::KEEP, comm);
}
diff --git a/src/Track/TrackRun.cpp b/src/Track/TrackRun.cpp
index 3be62f732d9aea237364b5e3dd15b1fefe0da1f0..9fe64b704eb7474ff28fdccd84721269de11d4ce 100644
--- a/src/Track/TrackRun.cpp
+++ b/src/Track/TrackRun.cpp
@@ -155,7 +155,7 @@ TrackRun* TrackRun::clone(const std::string& name) {
}
void TrackRun::execute() {
-
+
const int currentVersion = ((OPAL_VERSION_MAJOR * 100) + OPAL_VERSION_MINOR) * 100;
if (Options::version < currentVersion) {
@@ -241,13 +241,14 @@ void TrackRun::execute() {
macrocharge_m = beam->getChargePerParticle();
macromass_m = beam->getMassPerParticle();
- double Qtot = macrocharge_m * beam->getNumberOfParticles();
+ // double Qtot = macrocharge_m * beam->getNumberOfParticles();
/*
Here we can allocate the bunch
*/
-
+
+
// There's a change of units for particle mass that seems strange -> gives consistent Kinetic Energy
bunch_m = std::make_shared<bunch_type>(macrocharge_m,
beam->getMass()*1e9*Units::eV2MeV,
@@ -311,9 +312,13 @@ void TrackRun::execute() {
throw OpalException("Distribution::create", "Unknown \"TYPE\" of \"DISTRIBUTION\"");
}
+ *gmsg << "* About to create particles ..." << endl;
+
sampler_m->generateParticles(Np, nr);
- /*
+ *gmsg << "* Particle creation done" << endl;
+
+ /*
reset the fieldsolver with correct hr_m
based on the distribution
*/
@@ -361,8 +366,6 @@ void TrackRun::execute() {
Attributes::getBool(itsAttr[TRACKRUN::TRACKBACK]), Track::block->localTimeSteps,
Track::block->zstart, Track::block->zstop, Track::block->dT);
- *gmsg << "* Parallel Tracker created ... " << endl;
-
itsTracker_m->execute();
/*