From 53fa8ce22653adae83e04c1819274ac3ff5c070e Mon Sep 17 00:00:00 2001
From: Andreas Adelmann <andreas.adelmann@psi.ch>
Date: Sat, 22 Jun 2024 00:05:01 +0200
Subject: [PATCH] fix couplong constants
---
src/PartBunch/PartBunch.cpp | 71 ++++++++++++++++++++++++++++++++-----
src/PartBunch/PartBunch.hpp | 17 +++++++--
2 files changed, 78 insertions(+), 10 deletions(-)
diff --git a/src/PartBunch/PartBunch.cpp b/src/PartBunch/PartBunch.cpp
index 938fd2d8d..cf7b25ffd 100644
--- a/src/PartBunch/PartBunch.cpp
+++ b/src/PartBunch/PartBunch.cpp
@@ -2,6 +2,62 @@
#include <boost/numeric/ublas/io.hpp>
#include "Utilities/Util.h"
+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,10 +153,7 @@ void PartBunch<double,3>::calcBeamParameters() {
rmax_m(i) = rmax[i];
rmin_m(i) = rmin[i];
}
-
ippl::Comm->barrier();
-
-
}
template<>
@@ -285,9 +338,7 @@ void PartBunch<double,3>::computeSelfFields() {
//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),
@@ -297,18 +348,22 @@ void PartBunch<double,3>::computeSelfFields() {
},
Kokkos::Sum<double>(localDensity2) );
ippl::Comm->reduce(localDensity2, Density2, 1, std::plus<double>());
-
+
+ double Npoints = nr_m[0] * nr_m[1] * nr_m[2];
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);
}
diff --git a/src/PartBunch/PartBunch.hpp b/src/PartBunch/PartBunch.hpp
index 75a339f6b..ad4740ecf 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;
@@ -383,10 +391,11 @@ public:
*/
double getCouplingConstant() const {
- return 1.0;
+ return couplingConstant_m;
}
-
+
void setCouplingConstant(double c) {
+ couplingConstant_m = c;
}
void calcBeamParameters();
@@ -850,6 +859,10 @@ public:
double calculateAngle(double x, double y) {
return 0.0;
}
+
+ // Sanity check functions
+ void spaceChargeEFieldCheck();
+
};
/* \todo
--
GitLab