Code indexing in gitaly is broken and leads to code not being visible to the user. We work on the issue with highest priority.

Skip to content
Snippets Groups Projects
Commit 2ee03f85 authored by Sadr Mohsen's avatar Sadr Mohsen
Browse files

debugged. works on cpu

parent a4af3b20
No related branches found
No related tags found
1 merge request!15Multi variate normal
Hide whitespace changes
Inline Side-by-side
input-files/test-0004.in 0 → 100644
+ 59
− 0
View file @ 2ee03f85
/*
AWA Gun, configured in this input file for fast execution
(low particle number & low fieldsolver resolution, large
timestep and only one energy bin).
*/
OPTION, PSDUMPFREQ = 10; // 6d data written every 10th time step (h5).
OPTION, STATDUMPFREQ = 1; // Beam Stats written every time step (stat).
OPTION, BOUNDPDESTROYFQ = 10; // Delete lost particles, if any out of 10 \sigma
OPTION, AUTOPHASE = 4; // Autophase is on, and phase of max energy
// gain will be found automatically for cavities.
OPTION, VERSION=10900;
Title, string="AWA Photoinjector";
Value,{OPALVERSION};
// ----------------------------------------------------------------------------
// Global Parameters
REAL rf_freq = 1.3e9; // RF frequency. (Hz)
REAL n_particles = 1E6; // Number of particles in simulation.
REAL beam_bunch_charge = 1e-9; // Charge of bunch. (C)
//----------------------------------------------------------------------------
//Initial Momentum Calculation
REAL Edes = 1; //initial energy in GeV
REAL gamma = (Edes+EMASS)/EMASS;
REAL beta = sqrt(1-(1/gamma^2));
REAL P0 = gamma*beta*EMASS; //inital z momentum
//Printing initial energy and momentum to terminal output.
value, {Edes, P0};
D1: DRIFT, L = 1.0, ELEMEDGE = 0.0;
D2: DRIFT, L = 1.0, ELEMEDGE = 1.0;
D3: DRIFT, L = 1.0, ELEMEDGE = 2.0;
D4: DRIFT, L = 1.0, ELEMEDGE = 13.0;
L1: Line = (D1,D2);
L2: Line = (D3,D4);
BEAM1: BEAM, PARTICLE = ELECTRON, pc = P0, NPART = n_particles,
BFREQ = rf_freq, BCURRENT = beam_bunch_charge * rf_freq * 1E6, CHARGE = -1;
FS1: Fieldsolver, NX=8, NY=8, NZ=8, TYPE=FFT, PARFFTX = true, PARFFTY = true, PARFFTZ = false,
BCFFTX=OPEN, BCFFTY=OPEN, BCFFTZ=OPEN,
BBOXINCR = 1, GREENSF = INTEGRATED;
Dist1: DISTRIBUTION, TYPE= MULTIVARIATEGAUSS, SIGMAX=1e-1, SIGMAY=1e-2, SIGMAZ=1e-3, SIGMAPX=1e-2, SIGMAPY=1e-3, SIGMAPZ=2e-3, CORR = {1e-5,2e-5,0,0,0,0,0,0,0,0,0,0,0,0,0};
myLine: Line = (D1,D2,D3,D4);
TRACK, LINE = myLine, BEAM = BEAM1, MAXSTEPS = 1, DT = {1e-12}, ZSTOP=1.0;
RUN, METHOD = "PARALLEL", BEAM = BEAM1, FIELDSOLVER = FS1, DISTRIBUTION = Dist1;
ENDTRACK;
QUIT;
src/Distribution/Distribution.cpp
+ 14
− 20
View file @ 2ee03f85
......@@ -225,10 +225,7 @@ void Distribution::setDistParametersMultiVariateGauss() {
cutoffR_m = 3.;
cutoffP_m = 3.;
setSigmaR_m();
setSigmaP_m();
std::cout << "\n";
// initialize the covariance matrix to identity
for (unsigned int i = 0; i < 6; ++ i) {
for (unsigned int j = 0; j < 6; ++ j) {
if (i==j)
......@@ -243,8 +240,8 @@ void Distribution::setDistParametersMultiVariateGauss() {
setSigmaP_m();
for (unsigned int i = 0; i < 3; ++ i){
correlationMatrix_m[2*i ][2*i ] = sigmaR_m[i];
correlationMatrix_m[2*i+1][2*i+1] = sigmaP_m[i];
correlationMatrix_m[2*i ][2*i ] = sigmaR_m[i]*sigmaR_m[i];
correlationMatrix_m[2*i+1][2*i+1] = sigmaP_m[i]*sigmaP_m[i];
}
std::vector<double> cr = Attributes::getRealArray(itsAttr[DISTRIBUTION::CORR]);
......@@ -266,20 +263,8 @@ void Distribution::setDistParametersMultiVariateGauss() {
"Inconsistent set of correlations specified, check manual");
}
}
else{
}
avrgpz_m = 0.0;
std::cout << "\n";
for (unsigned int i = 0; i < 6; ++ i) {
for (unsigned int j = 0; j < 6; ++ j) {
std::cout << " " << correlationMatrix_m[i][j];
}
std::cout << "\n";
}
}
void Distribution::createDistributionGauss(size_t numberOfParticles, double massIneV, ippl::ParticleAttrib<ippl::Vector<double, 3>>& R, ippl::ParticleAttrib<ippl::Vector<double, 3>>& P, std::shared_ptr<ParticleContainer_t> &pc, std::shared_ptr<FieldContainer_t> &fc, Vector_t<double, 3> nr) {
......@@ -309,11 +294,20 @@ void Distribution::printDistMultiVariateGauss(Inform& os) const {
os << "* SIGMAPX = " << sigmaP_m[0] << " [Beta Gamma]" << endl;
os << "* SIGMAPY = " << sigmaP_m[1] << " [Beta Gamma]" << endl;
os << "* SIGMAPZ = " << sigmaP_m[2] << " [Beta Gamma]" << endl;
os << "* input cov matrix = ";
for (unsigned int i = 0; i < 6; ++ i) {
for (unsigned int j = 0; j < 6; ++ j) {
os << correlationMatrix_m[i][j] << " ";
}
os << endl << " ";
}
}
void Distribution::setAttributes() {
setSigmaR_m();
setSigmaP_m();
// setSigmaR_m();
// setSigmaP_m();
setDist();
}
void Distribution::setDist() {
......
src/Distribution/MultiVariateGaussian.hpp
+ 20
− 18
View file @ 2ee03f85
......@@ -48,16 +48,13 @@ public:
for (unsigned int i = 0; i < 6; i++) {
for (unsigned int j = 0; j <= i; j++) {
sum = 0.0;
for (unsigned int k = 0; k < j; k++) {
sum += L_m[i][k] * L_m[j][k];
}
if (j == i) {
for (unsigned int k = 0; k < j; k++) {
sum += L_m[j][k] * L_m[j][k];
}
L_m[j][j] = Kokkos::sqrt(cov_m[j][j] - sum);
}
else {
for (unsigned int k = 0; k < j; k++) {
sum += L_m[i][k] * L_m[j][k];
}
L_m[i][j] = (cov_m[i][j] - sum) / L_m[j][j];
}
}
......@@ -97,10 +94,15 @@ public:
}
for(int i=0; i<3; i++){
normRmin_m(i) = min_m(2*i);
normRmax_m(i) = max_m(2*i);
normPmin_m(i) = min_m(2*i+1);
normPmax_m(i) = max_m(2*i+1);
normRmin_m(i) = min_m(2*i)/opalDist_m->getSigmaR()[i];
normRmax_m(i) = max_m(2*i)/opalDist_m->getSigmaR()[i];
normPmin_m(i) = min_m(2*i+1)/opalDist_m->getSigmaP()[i];
normPmax_m(i) = max_m(2*i+1)/opalDist_m->getSigmaP()[i];
rmin_m(i) /= opalDist_m->getSigmaR()[i];
rmax_m(i) /= opalDist_m->getSigmaR()[i];
pmin_m(i) /= opalDist_m->getSigmaP()[i];
pmax_m(i) /= opalDist_m->getSigmaP()[i];
}
}
......@@ -163,26 +165,26 @@ public:
samplingP.generate(Pview, rand_pool64);
Matrix_t L;
for (unsigned int i = 0; i < 6; i++) {
for (unsigned int j = 0; j < 6; j++) {
for (unsigned i = 0; i < 6; ++i){
for (unsigned j = 0; j < 6; ++j){
L[i][j] = L_m[i][j];
}
}
}
// correlate samples by multiplying them into L
Kokkos::parallel_for( nlocal,KOKKOS_LAMBDA(const int k) {
double vec_old[6], vec[6];
double vec_old[6], vec[6] = {0.0};
for (unsigned i = 0; i < 3; ++i){
vec_old[2*i] = Rview(k)[i];
vec_old[2*i+1] = Pview(k)[i];
vec[2*i] = 0.0;
vec[2*i+1] = 0.0;
}
for (unsigned i = 0; i < 6; ++i){
for (unsigned j = 0; j < 6; ++j){
for (unsigned j = 0; j < i+1; ++j){
vec[i] += L[i][j]*vec_old[j];
}
}
for (unsigned i = 0; i < 3; ++i){
Rview(k)[i] = vec[2*i];
Pview(k)[i] = vec[2*i+1];
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment